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Abstract:  Results  from  the  efforts  in  meshfree  method  development  for 
fragment-impact  modeling  during  April  2009  to  March  2011  are  described 
in  this  report.  These  efforts  focused  on  an  enhanced  semi-Lagrangian 
reproducing  kernel  particle  method  formulation  for  modeling  large  material 
deformation  and  damage  mechanisms,  multiscale  homogenization  based  on 
the  energy  bridging  theory  pertinent  to  fragment  penetration  modeling  of 
concrete  materials,  and  enhanced  kernel  contact  algorithms  to  model  multi¬ 
body  contact  applicable  to  penetration  problems.  Several  benchmark 
problems  associated  with  contact-impact  simulations  as  well  as  multiscale 
modeling  of  material  damage  have  been  performed  to  examine  the 
effectiveness  of  the  developed  computational  methods.  These  newly 
developed  computational  formulations  and  the  associated  numerical 
algorithms  have  been  implemented  into  the  parallel  Nonlinear  Meshfree 
Analysis  Program  (NMAP)  code.  The  numerical  formulations  and  NMAP 
code  implementation  accomplished  under  this  effort  have  been  evaluated 
through  18  verification  and  validation  (V&V)  tests  based  on  penetration 
experiments  conducted  at  The  U.S.  Army  Engineer  Research  and 
Development  Center  (ERDC).  The  latest  version  of  the  NMAP  code  has  been 
delivered  to  ERDC.  This  report,  together  with  the  2009  technical  report 
(Chen  et  al.  2009),  offers  the  theoretical  foundation  of  NMAP  for  general 
users. 
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1  Introduction 

Impact  and  penetration  processes  in  brittle  materials  involve  many  complex 
phenomena  such  as  multi-body  contact,  fast-evolving  weak  and  strong 
discontinuities  resulting  from  material  softening  and  failure,  and  extensive 
material  deformation  and  flow.  Lagrangian  mesh-based  methods  are 
challenged  in  these  problems  as  a  result  of  excessive  mesh  distortion  and 
entanglement.  Mesh  sensitivity  associated  with  material  instability  in  strain 
localization  and  shear  band  formation  is  an  additional  issue.  In  terms  of  the 
material  description,  the  common  phenomenological  approach  to 
continuum  damage  evolution  lacks  a  direct  link  to  the  driving  micro¬ 
structure  failure  processes.  This  results  in  reduced  accuracy  in  terms  of 
describing  material  failure  and  a  loss  of  uniqueness  in  strain  localization 
problems. 

To  improve  the  capability  of  modeling  penetration  in  brittle  geomaterials, 
a  multiscale  reproducing  kernel  particle  method  (RKPM)  (Chen  et  al. 

1996)  formulation  and  a  new  damage  material  model  for  impact  and 
penetration  were  developed.  Major  project  developments  have  included: 

•  Enhanced  semi-Lagrangian  reproducing  kernel  (RK)  formulation  for 
modeling  excessive  material  deformation  and  damage  (Guan  et  al.,  in 
preparation) 

•  Multiscale  homogenization  based  on  the  energy  bridging  theory 
pertinent  to  fragment  penetration  modeling  of  concrete  materials  (Ren 
et  al.,  in  preparation) 

•  Enhanced  kernel  contact  algorithms  to  model  multi-body  contact 
applicable  to  penetration  problems 

•  Implementation  of  new  formulations  into  Nonlinear  Meshfree  Analysis 
Program  (NMAP)  parallel  code  (Chi  et  al.,  in  preparation) 

Organization  of  the  report  is  as  follows.  The  semi-Lagrangian  RK 
formulation  for  impact  and  penetration  is  presented  in  Chapter  2.  The 
multiscale  formulation  for  microcrack-informed  damage  evolution  is 
presented  in  Chapter  3.  Recent  enhancements  of  the  kernel  contact 
algorithm  are  given  in  Chapter  4.  Results  from  a  V&V  study  on  penetration 
modeling  with  the  NMAP  code  are  given  in  Chapter  5.  Conclusions  are 
given  in  Chapter  6. 
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2  Semi-Lagrangian  RKPM  for  Contact- 
Impact  Modeling 

For  large  deformation  problems  involving  severe  material  damage  and 
separation,  the  deformation  gradient  is  no  longer  positive  definite  at  all 
material  points  and  the  Total  Lagrangian  formulation  is  no  longer 
applicable.  Therefore,  to  address  the  limitation  of  Lagrangian  discretization 
for  modeling  materials  with  loss  of  positive  definiteness,  a  semi-Lagrangian 
RKPM  based  on  an  updated  Lagrangian  formulation  has  been  proposed 
(Guan  et  al.  2009;  Guan  et  al.,  in  preparation).  In  semi-Lagrangian  RKPM, 
the  mathematical  points  are  attached  to  the  material  points  (and  thus, 
Lagrangian),  while  the  supports  of  reproducing  kernel  shape  functions  do 
not  necessarily  deform  with  the  materials  (and  thus  semi-Lagrangian).  In 
this  approach  the  neighbor  points  are  redefined  during  the  deformation 
process  to  account  for  large  material  flow  and  formation  of  free  surfaces, 
and  the  approximation  is  constructed  at  the  current  configuration.  In  this 
work,  the  semi-Lagrangian  RKPM  framework  is  introduced  for  contact- 
impact  modeling.  The  stability  conditions  of  this  approach  for  explicit  time 
integration  are  also  presented. 

Reproducing  kernel  approximation 

Two  most  commonly  used  approximation  methods  in  meshfree  formula¬ 
tions  are  the  moving  least-squares  (MLS)  approximation  and  the  RK 
approximation.  Different  from  the  finite  element  method  (FEM)  where  the 
approximation  functions  are  constructed  in  the  mesh  mapped  to  a  natural 
coordinate,  MLS  and  RK  approximations  are  constructed  based  on  nodal 
positions  in  the  Cartesian  coordinate.  This  chapter  gives  a  review  of  the  RK 
approximation,  which  is  the  foundation  of  RKPM. 

Consider  a  closed  domain  Q  =  Q  u  dQ. ,  where  Q  is  the  open  domain  and 
3Q  is  the  boundary  of  Q ,  discretized  by  a  set  of  nodes  {xj ,  x2 ,  •  •  • ,  xNP } , 

x,  e  Q  ,  /  =  1, 2, . . NP  and  NP  is  the  number  of  points.  The  RK  approxi¬ 
mation  of  a  function  u,  denoted  by  uh ,  is  expressed  as  (Chen  et  al.  1996)  in 
Equation  1. 


NP 

1=1 


(1) 
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where  T7  (x)  is  the  RK  shape  function,  and  dl  is  the  corresponding 
coefficient.  The  RK  shape  function  is  constructed  with  the  following  form 

xPI(x)  =  C(x;x-xI)(pa(x-xI)  (2) 

where  x7  is  the  nodal  position  vector,  <pa  (x  -  x7 )  is  the  kernel  function,  and 
C  (x;  x  -  x, )  is  the  correction  function. 

The  kernel  function  cpa[x-xI)  is  a  compactly  supported  positive  function 

|<pa(x-x/)>°,  ||x  Xj- 1  /  a  <  1 

[<pa  (x  —  Xj )  =  0,  ||  x  —  Xj  ||  /  a  >  1 

where  a  is  the  measure  of  support  of  <pa  (x  -  x7 ) .  The  kernel  function 

expressed  in  Equation  3  has  a  spherical  support  with  radius  a.  A  kernel 
function  with  a  rectangular  or  cubic  support  can  be  constructed  by 
multiplication  of  one-dimensional  (l-D)  kernel  functions  (Chen  et  al. 

1996).  The  kernel  function  cpa  (x  -  x7 )  defines  the  locality  and  order  of 

continuity  of  the  approximation.  The  correction  function  C  (x;  x  -  x7 )  is 
the  combination  of  complete  nth  order  monomials 

C(x;x-xi)=  (*)(*1  —  */l)  (*2  —  *12)  (*3  —  *73) 

i+j+k= 0  (4) 

=  tfr(x-x7)6(x) 


HT  (x-Xj)  = 


X1~XI1 


X2  X 1 2 


X3  X 1 3 


q-x71)2  ...  (x3-x73)nl  (5) 


where  bijk  (x)  are  the  coefficients  of  the  basis  functions,  and  />(x)  and 
H  (x  -  x7 )  are  vectors  of  the  coefficients  and  monomial  basis  functions, 
respectively.  The  coefficient  vector  b(x)  is  solved  by  enforcing  the  exact 
reproduction  of  the  monomial  bases  up  to  the  nth  order. 

NP 

(*Kixi2*j3  =  *1*2*3  i  +  j  +  k  =  0,1  ,...,n  (6) 

1=1 


Equation  6  can  be  transformed  to  the  following. 
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J2WI  ( X)(X1  -  Xl/  )’  (X2  -  X2,  )J  (X3 


1=1 


-X3,)A'  5i08j05k0 


i+j  +  k  =  0,1,..., /i  (7) 


By  substituting  Equation  2  and  4  into  Equation7,  the  coefficient  vector 
b(x)  is  obtained  as 


b(x)=  M~1(x)H{0) 


(8) 


where 


NP 

M{X)  =  I2H{X-Xj)HT{X-Xj)(Pa{X-Xj)  /Q. 

J=1  w) 

ht(o)  =  [ 1  0  0  0  0  ...  0] 

Finally,  the  RK  shape  function  is  obtained  as 

Vi(jc)=  HT  (0)M~1(x)H(x-xI)q)a(x-xI)  (10) 

Figure  1  shows  the  contour  plot  of  a  two-dimensional  (2-D)  RK  shape 
function  with  rectangular  support.  The  kernel  support  is  shown  on  the  left, 
thus  defining  the  domain  of  influence,  and  the  shape  function  is  shown  on 
the  right,  constructed  using  a  cubic  B-spline  kernel  function  and  linear 
bases.  Figure  2  shows  the  comparison  of  RKPM  discretization  with  circular 
support  and  an  FEM  triangular  mesh  using  the  same  set  of  points.  The 
domain  of  influence  of  each  FEM  node  is  determined  by  the  neighboring 
connected  elements,  whereas  the  domain  of  influence  of  the  RK  shape 
function  is  defined  by  the  support  of  the  kernel  function.  While  in  RKPM 
discretization,  some  domains  of  influence  are  extended  outside  of  the 
physical  boundary,  the  reproducing  conditions  enforced  in  Equation  6 
guarantee  the  order  of  accuracy  for  all  veQ.  This  extended  boundary  layer 
in  RKPM  needs  to  be  considered  in  contact  problems.  However,  it  serves  as 
an  “insulation  layer”  to  ensure  impenetration  conditions  in  the  normal  con¬ 
tact  similar  to  the  function  of  a  “gap  element”  in  the  finite  element  setting. 
More  details  will  be  discussed  in  the  kernel  contact  algorithms  in  Chapter  4. 
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Figure  1.  Contour  for  2-D  RKPM  shape  function  with  rectangular  support: 
(a)  domain  of  influence,  (b)  RK  shape  function. 


(a) 


(b) 


Figure  2.  Comparison  of  FEM  and  RKPM  discretization  and  domains  of  influence: 

(a)  FEM  discretization  and  (b)  RKPM  discretization.  (The  domain  of  influence 
of  one  node  is  marked  in  grey  color  as  an  example.) 

Semi-Lagrangian  reproducing  kernel  discretization 

Governing  equations 

For  modeling  of  fragment-impact  processes,  a  semi-Lagrangian  RK 
discretization  is  introduced  to  the  equation  of  motion.  Start  with  an 
updated  Lagrangian  formulation  where  the  current  configuration  is  the 
referenced  configuration,  and  introduce  a  semi-Lagrangian  RK 
approximation  constructed  in  the  current  configuration  to  the  updated 
Lagrangian  variational  equation.  Let  X be  the  material  coordinate 
representing  the  initial  position  of  a  material  point,  and  x  be  the  current 
position  of  the  material  point  X  in  the  current  configuration  with  domain 
Qv ,  essential  boundary  dCl8x ,  and  natural  boundary  c)Q!'x  .  The  weak  form  of 

the  equation  of  motion  is: 
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L  +  L  Su^jdQ.  =  f  Su,b,da  +  f  Su,h,dr  (11) 


where: 

ul  =  displacement 

ry  =  Cauchy  stress,  u(iJ)  =  ( dui  /  dxj  +  du}  /  dxt)  /  2 

p  =  density 
bt  =  body  force 
ht  =  surface  traction 

In  the  pure  Lagrangian  RKPM  formulation,  the  Lagrangian  RK  shape 
functions,  VF/  (X) ,  are  constructed  using  the  material  coordinates  in  the 

initial  configuration.  The  discretization  of  Equation  11  by  the  Lagrangian 
RK  approximation  requires  taking  the  spatial  derivatives  of  the  Lagrangian 
RK  shape  function,  VF/  (X) ,  which  is  obtained  by  the  chain  rule  as 

d%UQ  d%pQdXj  dVjjX)  ! 
dxi  dXj  0xi  dXj  ji 

Here,  the  deformation  gradient  F  is  first  computed  using  the  material 
spatial  derivatives  of  the  Lagrangian  RK  shape  functions,  and  F~'  is  then 
obtained  by  taking  the  inversion  of  F  (instead  of  computing  F~l  directly). 
However,  the  Lagrangian  formulation  breaks  down  when  the  inverse  of  F  is 
no  longer  available.  This  may  occur,  for  example,  when  large  deformation 
leads  to  a  non-positive  definite  F  or  material  separation  takes  place.  Thus, 
in  this  research,  the  semi- Lagrangian  RKPM  formulation  (Guan  et  al.  2009) 
used  for  fragment-impact  problems. 

Semi-Lagrangian  discretization 

The  Lagrangian  formulation  breaks  down  when  mapping  between  the 
initial  and  current  configurations  is  no  longer  one-to-one.  This  occurs 
under  conditions  such  as  new  free  surface  formation  (i.e.,  material 
separation)  or  free  surface  closure,  which  commonly  exist  in  penetration 
processes.  Chen  and  Wu  (2007)  proposed  a  semi-Lagrangian  formulation 
to  overcome  the  shortcoming  of  the  Lagrangian  formulation. 
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In  the  semi-Lagrangian  formulation,  the  RKPM  points  follow  the  material 
motion,  while  the  distance  measure  z  =  x  -x(X  nt)  in  the  kernel  function 

is  defined  in  the  deformed  configuration.  Under  this  construction,  the 
kernel  support  of  the  semi-Lagrangian  kernel  function  does  not  deform 
with  the  material  motion.  The  comparison  of  the  supports  of  Lagrangian 
and  semi-Lagrangian  kernels  at  undeformed  and  deformed  states  is  shown 
in  Figure  3.  The  Lagrangian  kernel  supports  cover  the  same  group  of 
material  nodes  before  and  after  deformation,  while  the  semi-Lagrangian 
kernel  supports  cover  a  different  group  of  nodes  after  deformation. 


Figure  3.  Comparison  of  Lagrangian  and  semi-Lagrangian  kernel  supports  in  undeformed  and 
deformed  configurations:  (a)  undeformed  configuration,  (b)  Lagrangian  kernel  deformed  with 
the  material  in  the  deformed  configuration,  and  (c)  semi-Lagrangian  kernel  in  the  deformed 

configuration. 


The  semi-Lagrangian  RK  shape  function  formulated  in  the  current 
configuration  is  expressed  as 

^I(x)  =  HT{x-x(XI,t))b(x)cpa{x-x(XI,t))  (13) 


The  coefficient  vector  b(x)  is  solved  by  imposing  the  semi-Lagrangian 
reproducing  conditions. 


NP 

YJ'PI{x)x1(XI,t)1  x2(XI,t)Jx3(XI,tf  =x[xj2x\,  i  +  j  +  k  =  0,l,-,n  (14) 
1=1 

Solving  6(x)  from  Equation  14  and  substituting  it  into  Equation  13  yields 
the  following  semi-Lagrangian  RK  shape  function 


(*)  =  C(x;x-x(Xn  t))(pa  (x  -  x(Xj  ,f)) 


(15) 
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where 

C(x\x  -  x  (XItt))  =  HT  (t 0)M -1  (. x)H(x  -x(Xnt))  (16) 

NP 

M(x)  =  Y/H(x-x(XI,t))HT(x-x(XI,t))(Pa(x-x(XI,t))  (17) 
1=1 

HT(x-x(Xnt))  =  \l  x1-x1(XI,t)  x2-x2(Xnt )  (x3-x3(X;,t))"]  (18) 

Let  the  velocity,  v. ,  be  approximated  by  semi-Lagrangian  RK  shape 
functions. 

NP 

v^x,t)  =  ^%(x)vu(t)  (19) 

1=1 

The  corresponding  semi-Lagrangian  approximation  of  acceleration  is 
given  as 

NP 

u-(x,t)  =  v-(x,t)  =  ^(^Wr^O  +  'fjWr^f))  (20) 

1=1 

where  T^x)  is  the  correction  due  to  the  time  rate  of  the  semi-Lagrangian 
kernel  <pa 

V;(x)  =  C(x;x-x(XI,t))<pa(x-x(XI,t))  (21) 


x  —  x(X,,o  ,  x  —  x(X,,t)  «-(v  —  vT 
<3 Pa  - "  =(Pa  - "  — - 

a  a  a 


•  d(  ) 

where  (  )  =  ’  is  the  material  time  derivative,  and 


X  —  X^Xjyt) 

x  —  x(Xj,t) 


(22) 


(23) 


Substituting  Equation  20  into  Equation  li  yields  the  following  semi- 
Lagrangian  discrete  equation 
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Mv  +  Nv  =  fext-fint  (24) 

where 

M,J=fp<r,(x)vJ(x)/dn  (25) 

h. 

Nu  =  ffiH',(x)'}C(x)lda  (26) 

n 

/,“  =  f  Bjsda  <27) 

% 

/“'  =  f  WITbdQ+  J  WjhdT  (28) 

nx  onx 


where  B,  is  the  gradient  matrix  associated  with  uU  j) ,  E  is  the  stress 
vector  associated  with  rjy ,  and  b  and  h  are  body  force  and  surface  traction 
vectors,  respectively. 

Stabilized  Non-conforming  Nodal  Integration  (SNNI) 

Domain  integration  in  Galerkin  meshfree  methods  requires  special 
attention  as  there  is  no  mesh  in  the  discretization.  Gauss  integration 
requires  a  background  grid  and  introduces  significant  integration  errors  if 
the  kernel  supports  do  not  match  with  the  integration  grid.  Nodal 
integration  with  stabilization,  such  as  the  stabilized  conformation  nodal 
integration  (SCNI)  (Chen  et  al.  2001;  Chen  et  al.  2002)  was  proposed  as  an 
alternative  to  the  Gauss  integration.  In  SCNI,  nodal  strain  smoothing  on  a 
conforming  smoothing  domain  surrounding  each  node  is  introduced  to 
achieve  stability  and  optimal  convergence.  The  strain  smoothing  at  the 
node  L  is  calculated  by 

«(/k)=-r  Jeedn=Trf{ui,i+uj.i)dn  <29> 

L  aL  nL 


where: 
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£U 


and 


components  of  the  regular  strain  and  smoothed 

strain  tensors,  respectively 
components  of  the  displacement  vector 

area  (or  volume)  of  the  conforming  smoothing 
domain  associated  with  node  L 


By  introducing  the  divergence  theorem,  the  previous  equation  is 
transformed  to  a  boundary  integral  to  yield 


-zrfK+u»)da 


=jrJ(u<ni+uin<)dr 

L  anL 

*  NP  NP 

=  JT  f  +EWJ  (30) 

dnL  1=1  1=1 

NP  _  _ 

=  E(Vhz+Vhy) 

i=i 


where 


*«  =  4- /*,(*)", dr 

L  anL 


(31) 


In  Equation  31,  by  is  the  smoothed  gradient  of  the  shape  function.  Strain 
smoothing  on  a  conforming  smoothing  domain,  that  is,  uQL  =Q,  is  the 

requirement  to  satisfy  the  integration  constraint  for  optimal  convergence 
(Chen  et  al.  2001).  One  choice  for  generation  of  the  conforming  smoothing 
domains  is  the  Voronoi  diagram,  as  shown  in  Figure  4.  However,  SCNI 
comes  with  the  cost  of  constructing  the  conforming  smoothing  domains, 
and  it  is  particularly  tedious  in  the  semi-Lagrangian  discretization  where 
smoothing  domain  reconstruction  at  every  time-step  is  needed.  For 
penetration  problems,  the  damage  evolution  and  the  associated  new 
surface  formation  further  complicate  the  Voronoi  cell  generation.  In  this 
work,  the  stabilized  non-conforming  nodal  integration  (SNNI)  (Chen  et  al. 
2006;  Chen  and  Wu  2007)  is  introduced.  In  this  approach,  the  conforming 
requirement  in  the  nodal  strain  smoothing  domain  is  not  enforced,  that  is, 
as  shown  in  Figure  5,  where  simple  strain  smoothing  domains 

are  used.  RKPM  with  SNNI  remains  stable,  with  accuracy  comparable  to 
that  with  SCNI  (Chen  et  al.  2006;  Chen  and  Wu  2007). 
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Figure  4.  Strain  smoothing  domains  for  SCNI. 


Boundary  correction  of  SNNI  domain  integration 

How  the  strain  smoothing  domains  of  SNNI  should  be  properly  corrected 
on  the  problem  boundaries  are  shown  here.  Consider  solving  the  elastic 
wave  propagation  through  a  three-dimensional  (3-D)  rod  using  semi- 
Lagrangian  RKPM  with  SNNI  domain  integration.  The  properties  of  the  rod 
are:  length  =  iom(0<x<10),  square  cross-sectional  area  ( 0  <  y,  z  <  0.02 ), 
Young’s  modulus  E  =  10000 ,  Poisson’s  ratio  v  =  0 ,  and  density  p- 1 .  Here 
the  Poisson’s  effect  is  purposely  removed  for  the  wave  to  propagate  only  in 
the  axial  direction,  while  using  a  3-D  geometry  to  magnify  the  boundary 
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effect.  The  domain  is  discretized  by  100  x  3  x  3  points  in  x,  y  and  z  directions. 
The  boundary  and  initial  conditions  are 


Boundary  conditions 


u,  (x  =  0,t)  =  0 
OyUj  on  all  other  surfaces 


initial  condition 


vt  (x,t  =  0)  =  Si±  sin 
u;  (jc,f  =  0)  =  0 


XJl 


2  L 


(32) 


The  uncorrected  smoothing  domains  for  SNNI  are  first  considered  as 
shown  in  Figure  6.  The  displacement  and  velocity  histories  at  the  middle 
point  x=5  are  plotted  in  Figure  7.  It  is  observed  that  the  error  starts  to 
accumulate  when  the  wave  reflects  from  the  boundary,  as  shown  in  Figure 
7(a).  This  is  attributed  to  the  fact  that  smoothing  domains  of  the  boundary 
nodes  are  extended  outside  of  the  problem  domain.  This  is  corrected  by 
shifting  the  boundary  smoothing  domains  inward,  as  shown  in  Figure  8. 
The  numerical  solution  with  correction  of  the  boundary  smoothing 
domains  agrees  well  with  the  analytical  solution,  as  shown  in  Figure  7b. 


Temporal  stability  of  semi-Lagrangian  RKPM 

This  section  discusses  stability  of  the  semi-Lagrangian  RKPM  formulation. 
Under  the  semi-Lagrangian  RKPM  formulation,  the  RK  shape  functions 
are  constructed  at  every  time-step,  which  allows  the  particles  to  transport 
in  and  out  of  the  support  of  each  RK  shape  function.  This  gives  the 
flexibility  to  define  material  contact  and  separation  at  the  discrete  level. 
Under  this  construction,  the  material  time  derivative  of  semi-Lagrangian 
discretized  variables  is  calculated  in  accordance  with  Equations  20 
through  24,  and  the  discrete  equation  of  motion  is  given  in  Equation  24. 
Considering  a  central  difference  temporal  discretization,  the  following  full 
discrete  equation  is  given. 


Using  the  von  Neumann  stability  analysis  for  a  l-D  wave  equation 
discretized  by  semi-Lagrangian  RKPM,  the  amplification  factor  can  be 
obtained  as  follows. 


Consider  a  wave  form  of  the  discrete  solution  d"  =  [Xf  ek{l^] ,  where  k  is  the 

wave  number  and  2  is  the  amplification  factor.  From  Equation  33  the  the 
following  equation  is  obtained  for  the  amplification  factor  (following  Guan 
et  al.  2009). 
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Figure  6.  Side  view  of  the  smoothing  domains  on  the  left  end  of  the 
rod.  Solid  thick  lines:  problem  boundary;  thin  solid  lines:  boundaries  of  smoothing  domains. 


Figure  7.  Wave  propagation  solved  by  semi-Lagrangian  RKPM  with  SNNI 
domain  integration  method:  (a)  without  boundary  SNNI  smoothing  zone 
correction,  (b)  with  boundary  SNNI  smoothing  zone  correction. 


Figure  8.  Correction  of  the  boundary  smoothing  domains. 
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M (< in+1  -2 dn+dn~1)  =  (Atf  l^fext  )n+1  -  (/int )"  -  Nvn  J 

\  +  2S1+^S2=0 
A 

At2 

S1  =-1  +  2 — -c2  A  +pIAtAxA1 
Ax 

S2=  1  -  2p,  AxAfW*  (Ax)  ^ 

'T  (Ax)  =  C(x,  +  Ax;  Ax)  <4  (Ax)^/+1~Uj^ 


(33) 


(34) 


Here  the  coefficients  /t,  and  +  are  related  to  the  method  of  domain 

integration,  and  are  given  in  Guan  et  al.  (in  preparation)  for  the  case 
where  SNNI  is  introduced.  The  contour  plot  of  the  amplification  factor  is 
shown  in  Figure  9  for  the  case  where  the  kernel  function  support  size  is 
chosen  to  be  1.5  times  the  nodal  distance  (R=i.s).  The  temporal  stability 
requires  the  amplification  factor  to  be  bounded  by  unity,  |l|  <  1 ,  as  shown 

in  Figure  9,  where  the  horizontal  axis  is  the  normalized  time-step  and  the 
vertical  axis  is  the  velocity  gradient  between  two  adjacent  nodes.  The 
additional  parameter  a  shown  in  the  figure  is  the  ratio  between  the  size  of 
the  strain  smoothing  domain  and  the  nodal  distance,  Ax .  The  results  show 
that  the  stability  of  semi-Lagrangian  RKPM  is  affected  by  the  velocity 
gradient.  This  needs  to  be  considered  in  the  kernel  contact  algorithms,  as 
will  be  discussed  in  Chapter  4. 


\X\  for  SNNI  with  a  =  1.5 


Ax 


Figure  9.  Contour  plot  of  amplification  factor  of  semi-Lagrangian  discrete  equation 

with  SNNI. 
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Implementation 

In  this  section,  the  general  coding  structure  of  the  transient  semi- 
Lagrangian  RKPM  formulation  is  introduced.  The  flow  chart  is  shown  in 
Figure  10.  In  semi-Lagrangian  RKPM,  the  shape  functions  are  constructed 
in  the  current  configuration,  which  calls  for  efficient  search  of  neighbors.  In 
this  work,  a  simple  bucket  search  scheme  is  introduced.  Since  RK  shape 
functions  do  not  have  the  kronecker  delta  property,  VP/  (jc, )  *  5U ,  it  requires 

calculation  of  displacements  using  shape  functions  and  generalized 
coordinates  (Chen  and  Wang  2000).  By  using  nodal  integration  in  the  semi- 
Lagrangian  RKPM,  all  state  and  field  variables  are  stored  at  nodal  points. 


Figure  10.  Flow  chart  of  dynamic  semi-Lagrangian  RKPM  code  structure. 


A  major  effort  in  semi-Lagrangian  RKPM  is  the  neighbor  search  at  each 
time-step.  Consider  a  RK  approximation  of  u(x)  at  xL,  which  requires 

the  search  of  all  the  nodes  whose  kernel  supports  cover  xL  as  shown  in 
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Figure  n.  This  group  of  neighbor  nodes  for  xL  is  denoted  as 
N L  =  {./ 1  <pa  (.v^  -  x, )  *  0  j .  A  two-level  bucket  search  algorithm  is  intro¬ 
duced.  In  the  first  step,  nodes  located  inside  a  brick  domain  bounded  by 
xL-h  and  xL+h ,  where  It  =  ( hx ,  hy ,  h_  j  with  ht  properly  selected  to  be  suf¬ 
ficiently  larger  than  the  kernel  support,  are  marked  as  candidate  neighbors 
depicted  as  open  dashed  circles  in  Figure  12.  Nodal  distances  between  the 
candidate  neighbors  and  xL  are  then  calculated  to  identify  NL  as  shown  in 

closed  circles  in  Figure  12.  The  bucket  of  candidate  neighbors  is  updated 
according  to  the  local  velocity  field  in  the  bucket.  As  shown  in  Figure  13,  the 
update  time  is  defined  as  the  time  required  for  the  node  L  to  move  more 
than  the  minimum  support  size  of  the  candidate  neighbors.  This  requires 
the  velocity  of  the  node  L  calculated  by 

vUJ  =  2  (35) 

l€NL 


Figure  11.  Neighbors  of  point  L  marked  as  black  dot. 


•  •  •  • 


Figure  12.  Bucket  search  of  candidate  neighbors  of 
point  L  (open  dash  circle  points). 
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3  Multiscale  Microcrack  Informed  Damage 
Model 

Damage  processes  in  brittle  solids  are  driven  by  the  distribution  of 
microcracks  and  their  evolution.  Explicit  modeling  of  microcracks  in  brittle 
solids  is  computationally  infeasible,  while  classical  damage  models  are  phe¬ 
nomenological  in  nature  with  a  missing  link  to  microscopic  properties.  A 
micromechanics-based  approach  was  introduced  (Nemat-Nasser  and  Hori 
1999)  to  obtain  homogenized  material  properties  of  cracked  solids.  How¬ 
ever,  they  are  limited  to  certain  idealized  crack  configurations  and  loading 
conditions.  Alternatively,  asymptotic  expansion  offers  a  rigorous  means  for 
relating  physical  variables  defined  at  different  scales  (Bakhvalov  and 
Panasenko  1989;  Benssousan  1978;  Fish  1997;  Guedes  and  Kikuchi  1990). 
However,  the  key  step  in  an  asymptotic-type  method  is  solution  of  the 
characteristic  functions  in  the  microscopic  cell,  which  is  typically  very  time- 
consuming. 

The  present  work  aims  at  constructing  damage  models  based  on 
thermodynamics  of  “cracked”  microscopic  cells  and  the  corresponding 
“damaged”  macroscopic  continua.  Crack  evolution  in  the  microscopic  cell 
is  first  modeled  numerically.  Then,  corresponding  damage  evolution 
functions  in  the  continuum  are  constructed  through  a  Helmholtz  free 
energy  relationship  between  damaged  and  undamaged  homogenized 
continua.  This  approach  avoids  the  tedious  solution  of  characteristic 
functions  in  the  conventional  asymptotic  type-methods. 

Model  problem  and  homogenization  operators 

Start  with  a  two-scale  representation  of  the  model  problem.  A  hetero¬ 
geneous  solid  with  domain  Q  and  boundary  Y  containing  a  distribution  of 
microcracks  is  considered,  as  shown  in  Figure  14.  The  solid  is  subjected  to 
surface  traction  t  on  T,  and  prescribed  displacement  u  on  Tu  ,  Tr  UTh  =  T, 

and  is  assumed  to  undergo  static  elastic  deformation  without  body  force. 
For  a  given  material  point  in  the  macroscopic  solid,  it  corresponds  to  a 
microscopic  cell  microstructure  with  domain  Qv  containing  microcracks 

with  surface,  Yc .  Use  x  as  the  macroscopic  coordinate  and  y  as  the 
microscopic  coordinate. 
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Figure  14.  Microscopic  and  macroscopic  structures. 
The  boundary-value  problem  is  given  as 


V  •  o£  =  0  mfl 

(36) 

with  the  corresponding  boundary  conditions 

o£  •  n  =  t  on  Tt 

(37) 

u£  =  U  on  r„ 

(38) 

and  the  effective  traction  acting  on  the  crack  surface 

a£n  =  h  on  Tc  (39) 


where  c  is  the  total  stress  field,  uf  is  the  total  displacement  field, 
superscript  “s  ”  denotes  that  the  coarse  and  fine  scale  responses  are 
embedded  in  the  total  solution,  n  is  the  unit  outward  normal  vector  on  the 
boundary,  and  the  traction,  h ,  is  applied  to  the  union  of  crack  surfaces,  Yc 

Consider  a  linear  elastic  material  law 


oe  =  Ce :  ee 


(40) 
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where  CE  is  the  heterogeneous  constitutive  tensor,  and  e£  is  the  total 
strain  field,  where 


e£=|(v®u£+u£®v)  (41) 

Direct  simulation  of  the  total  scale  governing  Equations  36  to  39  is  time 
consuming  due  to  the  fine-scale  microscopic  defects  and  heterogeneities, 
and  an  attempt  is  made  to  conduct  homogenization  as  shown  in  Figure  15, 
where  the  homogenized  stress  and  strain  are  defined  by  the  microscopic 
cell. 


Figure  15.  Homogenization  of  a  microscopic  cell  with  fluctuating  fields. 


The  tractions  and  displacements  prescribed  on  the  outer  boundary  of  the 
microscopic  cell,  dQy ,  are  related  to  the  homogenized  stress  and  strain  in 

the  continuum  as 


a 


frf  (t"®x)ds 

y  d£ly 


(42) 


y  dn„ 


)  n  +  n  ( 


>u£)ds 


(43) 


where  Vy  =  |  <r/Q  is  the  volume  of  the  microscopic  cell.  Alternatively,  the 
averaged  stress  and  strain  in  the  microscopic  cell  are  defined  as 
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(o£\  =  —  f  o£dO 
X  1  Vy 

(es)  =  —  f  e£dO 
X  '  Vy 

where 


(44) 


(45) 


(46) 


The  following  equation  is  used  to  obtain  the  relationship  between 
homogenized  and  averaged  stresses  for  the  cracked  microscopic  cell 

V-(o£®x)  =  V-o£®x  +  a£-(V®x)  =  a£  (47) 


By  substituting  Equation  47  into  the  averaged  stress  definition  in  Equation 
44,  the  following  equation  is  used 


±f(t‘0X)ds- 

y 


— J  V-(o£  <g>x)  dO 

Vv  ay 

±f(t‘0x)ds 

y  rc 


=  O 


V 


■(f)  (t £  ®x)ds 


y  1; 


(48) 


The  second  term  on  the  right  hand  side  of  Equation  48  vanishes  due  to 
equilibrium  of  the  cohesive  stresses  on  the  crack  surface.  Hence,  the 
homogenized  stress  equals  to  the  averaged  stress  even  in  the  case  of  a 
cracked  microscopic  cell,  that  is 


o  =  (o£)  (49) 

Substituting  Equation  41  into  averaged  strain  defined  in  Equation  45  and 
considering  the  divergence  theorem,  the  following  is  created 
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}  =  —  J  e£dO  =  —  J(v®u£+u£®v)dO 


V  J  2V 

y  y  a 


j)  (u£  0n  +  n0usJds — — ^(u£  <g>n  +  n<g>us)ds  (50) 

on 


y  r 


2V 


-^(u£  ®n  +  n®u£)ds 


y  r. 


Consequently  the  relationship  is  obtained  between  homogenized  strain 
and  averaged  strain  as 


e  = 


u£  ®n  +  n®u£ 


(51) 


Here  it  is  shown  that  the  homogenized  strain  consists  of  the  average  strain 
and  the  additional  strain  induced  by  the  displacement  jump  across  the 
crack  surfaces.  In  other  words,  macroscopic  strain  cannot  be  directly 
obtained  by  the  averaging  of  microscopic  strain  when  microcracks  exist. 

To  obtain  the  homogenized  material  tensor,  the  fourth-order  influence 
tensor  is  defined  as 


e£(x)  =  [l-A£(x)]:e  (52) 

where  I  is  the  fourth-order  identity  tensor.  Substituting  Equation  52  into 
the  stress-strain  relation  in  Equation  40,  the  following  is  obtained 

o£  =C£  :e£  =C£  :[l-A£(x)]:e  (53) 

Using  the  equivalence  between  the  homogenized  stress  and  averaged 
stress,  obtained  is 


d  =  (c£  :[l-A£(x)]^:e  =  C:e  (54) 

From  the  previous  equation,  the  homogenized  material  response  tensor  is 
expressed  as 


C  =  (I-D):(C£) 


(55) 
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where 


D  =  (c':A-):(c')_1 


(56) 


Here  D  is  a  degradation  (damage)  tensor  expressed  by  the  influence 
tensor,  A£(x) . 

Remark  3.1 

The  fourth-order  influence  tensor  A*(x)  represents  the  interscale  relation 
between  properties  of  the  cracked  microstructures  and  properties  of  the 
homogenized  continua.  Micromechanics  offer  an  analytical  basis  for 
obtaining  the  influence  tensor,  for  example,  Mori-Tanaka  method  (Mori 
and  Tanaka  1973),  self-consistent  method  (Budiansky  and  O’Connell 
1976),  and  differential  scheme  (Norris  1985).  On  the  other  hand, 
asymptotic  expansion  based  homogenization  offers  a  general  approach  for 
calculating  the  influence  tensor,  which  utilizes  the  numerical  solution  of 
characteristic  functions  for  scale-coupling  in  the  microscopic  cell. 

Remark  3.2 

Using  the  asymptotic  expansion  based  method,  the  influence  tensor  can  be 
obtained  from  the  characteristic  tensor  %(y)  as  (see  Appendix  A  for 

detailed  derivations) 


A£  =  -VsyX(y)  (57) 

where  the  characteristic  tensor  /(y)  is  solved  from  the  following 
boundary-value  problem  defined  in  the  microscopic  cell. 

vy-(c£:v;x(y))  =  0  in  ny  (58) 

[c£ :  Vyx(y):  V£v[0](x)]-n  =  -[ce :  V£v[0](x)]  n  +  h  on  Tc  (59) 

Clearly  the  numerical  solution  of  the  previous  boundary- value  problem  for 
obtaining  the  third-order  tensor  /(y)  is  computationally  intensive.  This 

issue  is  addressed  by  introducing  the  energy  bridging  method  as  discussed 
in  the  next  section. 
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Energy  bridging  between  scales 

To  establish  the  relationship  between  the  microcracks  induced  material 
degradation  and  the  damaged  continuum  described  by  conventional 
damage  mechanics,  the  Helmholtz  free  energy  (HFE)  is  employed. 

m  1-  - 

W  =  -o\e  (60) 


According  to  the  second  principle  of  thermodynamics,  is 


dW 


de 


dw 

dD 


(61) 

(62) 


where  D  is  the  damage  tensor  and  Y  is  the  damage  energy  release  rate 
(DERR)  representing  a  driving  force  of  damage  evolution.  Consequently, 
the  damage  evolution  equation  is  expressed  by  the  evolution  potential,  ®  , 
as 


• 

D  =  A—  (63) 

oY 

where  X  is  the  consistency  parameter.  In  conventional  continuum  damage 
mechanics,  HFE  is  determined  experientially.  In  the  present  approach,  the 
HFE  is  obtained  through  homogenization  of  the  cracked  microstructure. 
The  microscopic  free  energy  is  defined  as 

'Fe=|o£:e£  (64) 


where  <f  and  e£'  are  the  microscopic  stress  and  strain.  Integrating 
microscopic  HFE  in  the  microscopic  cell  yields 
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f  Wedn  =  -  J ae :  e£  dQ. 

n y 

=  -  J-[v-(u£-a£)  +  V-(a£-u£)-u£-V-a£-V-a£-u£]cttl 

2  n9  2 

=  lfv(u--o-)da 


(65) 


Note  that  equilibrium  was  used  without  body  force.  Introducing  the 
divergence  theorem  in  Equation  65  results  in 


J  Wed£l  =  -  J  V-(u  eoe)dn 


1 

2 

1 

2 


(j)  ue  oe  nds-—  J u£  o£  n  ds 

dny  2  rc 

j)  u£  t£ds-^^u£  hds 


(66) 


Here  a  strain  driven  homogenization  is  considered,  where  the  displacement 
boundary  conditions  obtained  from  macroscopic  strain  are  applied  to  the 
microscopic  cell  to  compute  the  local  stress  field.  The  local  stress  field  is 
then  averaged,  which  is  the  same  as  the  homogenized  stress  a  according  to 
Equation49,  and  is  passed  back  to  the  coarse  scale.  The  prescribed 
boundary  displacements  on  the  microscopic  cell  are  obtained  from  the 
macroscopic  strain  by 


u£=e-x  on  d£2y  (67) 

where  dQy  is  the  outer  boundary  of  the  microscopic  cell.  Substituting 
Equation  67  into  the  first  term  on  the  right  hand  side  of  Equation  66  is 
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—  j)  u£  •  t£  ds  =  —  (e  •  x)  •  t£  ds 

2  dny  2  dny 


1 

2 


j)  (t£  ®x)  ds 

any 


e 


Vu 

—  —d :  e  =  F  ? 
2  y 


(68) 


Combining  Equations  66  and  68  is 


(69) 


The  right  hand  side  of  Equation  69  is  the  averaged  free  energy  of  the 
microscopic  cell,  while  the  left  hand  side  of  Equation  69  is  the  energy 
density  of  the  homogenized  material.  Thus,  Equation  69  concludes  that 
the  averaged  energy  of  the  microscopic  cell  is  equal  to  the  energy  density 
of  the  homogenized  material. 

Remark  3.3 

The  result  in  Equation  69  is  an  extension  of  Hill’s  theorem  (Hill  1963), 
with  addition  of  the  second  term  on  the  right  hand  side  resulting  from  the 
crack  discontinuous  displacement  field.  This  result  is  consistent  with  the 
one  obtained  by  (Belytschko  et  al.  2008)  where  the  power  expression  of 
the  energy  bridging  equation  was  introduced.  Further,  the  energy 
equivalence  in  Equation  69  can  also  be  obtained  by  the  asymptotic 
expansion  approach  as  detailed  in  Appendix  B. 

Characterization  of  parameters  in  continuum  damage  mechanics  by 
energy  bridging 

The  Helmholtz  free  energy  relationship  is  used  between  the  homogenized 
continuum  and  the  cracked  microstructure  to  derive  damage  parameters 
for  several  commonly  used  damage  models. 

One-parameter  scalar  damage  model 

The  one-parameter  damage  model  (Mazars  1984)  is  expressed  as 
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W  =  (l  -d)W0  (70) 

where  'k  is  Helmholtz  free  energy  computed  from  the  microscopic  cell  in 
Equation69,  and  To  is  the  effective  Helmholtz  free  energy  obtained  by 

To=|e:C0:e  (71) 

where  C0  is  the  homogenized  undamaged  material  tensor.  Then  the 
damage  scalar  is  obtained  by 


d  =  1 


W 

W 


(72) 


Two-parameter  damage  model 

The  two-parameter  damage  model  was  extensively  studied  and  widely 
adopted  for  description  of  more  complicated  damage  mechanisms  (Faria 
et  al.  1998;  Li  and  Ren  2009;  Wu  et  al.  2006).  For  geological  materials 
such  as  sand  or  soil,  the  volumetric-deviatoric  splitting  approach  was 
adopted  to  describe  damage  mechanisms  driven  by  pressure  and  shear.  On 
the  other  hand,  the  hydrostatic  tension-compression  decomposition 
approach  was  used  for  quasi-brittle  materials  such  as  concrete  and  rock 
(Faria  et  al.  1998),  which  is  the  case  demonstrated  herein.  The  framework 
of  the  two-parameter  damage  model  is  described  in  Appendix  C. 

The  initially  undamaged  macroscopic  stress  is  obtained  as 

d0  =  C0 :  e  (73) 

The  macroscopic  stress  can  be  decomposed  as 

d0  =  0+  +b0  (74) 

where  and  are  the  tensile  and  compressive  stresses,  respectively, 
defined  as 


(75) 
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o"=o0-o+  (76) 

where  <j.  and  p;  are  the  ith  eigenvalue  and  the  corresponding  eigenvector 
of  o0 ,  andH(  )is  the  Heaviside  function. 

The  macroscopic  Helmholtz  free  energy  is  then  expressed  as 

W  =  (l~ d+)X  +  (1  - dW0  (77) 

where  d+  and  d~  are  tensile  and  compressive  damage  parameters, 
respectively,  and  the  corresponding  expressions  of  the  effective  Helmholtz 
free  energy  are 


1  _ 

in*  —  •  p1  .  n± 

—  2  °°  '  *  °° 


(78) 


Considering  the  Clausius-Duhem  inequality  of  irreversible  thermo¬ 
dynamics,  the  following  can  be  derived  from  the  Helmholtz  free  energy  in 
Equation  77 


Y±  =  Vo 


(79) 


and 


d± 


dW  ,  AW 

()Y  ~  "AT* 


(80) 


where  T  is  the  Helmholtz  free  energy  calculated  from  the  microscopic  cell 
in  Equation  69.  By  using  the  finite  difference  operation  in  Equation  80, 
the  evolution  of  damage  parameters  is  obtained. 

Fully  tensorial  damage  model 

Derivation  of  the  tensorial  damage  tensor  based  on  the  influence  tensor 
Af(x)  according  to  Equation  56  is  limited  to  special  cases  if  A"'(x)  is 
obtained  by  the  micromechanical  methods,  and  is  costly  if  A£(x)  is 
calculated  using  the  asymptotic  expansion  method.  By  introducing 
bridging  based  on  Helmholtz  free  energy,  a  more  efficient  approach  for 
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obtaining  the  damage  tensor  is  introduced.  Recalling  the  definition  of  the 
macroscopic  Helmholtz  free  energy  in  Equation  60  the  following  appears 


W  =  |e:(l-D):C0:e 


(81) 


Performing  partial  differentiation  of  the  Helmholtz  free  energy  in 
Equation  81  with  respect  to  the  fourth-order  damage  tensor,  D ,  the 
damage  energy  release  rate  is  obtained  as 

dW  1  - 

Y  =  —  =  — e :  C0 :  e  (82) 

<9D  2  0 


By  taking  the  derivative  of  the  Helmholtz  free  energy  of  the  microscopic 
cell  in  Equation  69  with  respect  to  the  damage  energy  release  rate,  Y ,  the 
fourth-order  damage  tensor  is  obtained. 


dW  ,  AW 

_  I _ 

d\  AY 


(83) 


Here  a  finite  difference  approach  could  be  used  in  Equation  83  to  obtain 
the  fourth-order  damage  tensor  using  A4*  and  AY . 

Microscopic  cell  analysis 

The  essential  step  in  obtaining  the  damage  evolution  functions  in  the 
proposed  method  is  the  microscopic  cell  analysis.  In  the  present  study,  the 
microscopic  cell  is  made  of  an  isotropic  linear  elastic  material  with  a 
center  crack  subjected  to  boundary  conditions  as  shown  in  Figure  14.  The 
boundary  conditions  on  the  microscopic  cell  are  defined  based  on  the  type 
of  damage  mechanisms  to  be  captured.  The  crack  is  modeled  by  a  cohesive 
crack  model  with  linear  degradation  as  shown  in  Figure  14.  The  crack 
propagation  direction  is  determined  by  the  maximum  hoop  stress 
criterion. 


In  this  work,  we  consider  an  enriched  reproducing  kernel  particle  method 
(RKPM)  for  the  microscopic  cell  analysis.  The  microscopic  cell  solution  is 
then  used  in  the  homogenization  as  previously  described  for  characteriza¬ 
tion  of  the  damage  parameter  evolution. 
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In  RKPM  (Chen  et  al.  1996;  Liu  1995),  the  approximation  of  a  function  u  is 
expressed  as 


uh00  =  '%2fi(x)u1  (84) 

I 

where  4(x)  is  the  reproducing  kernel  (RK)  shape  function  and  n,  is  the 
corresponding  coefficient.  In  the  following  equations,  we  use  multi¬ 
dimensional  notation:  a  =  (a±,a2,- •  -ad) ,  |a|  =  y  ^  a- ,  xa  =  x"1  •  ■■xd' , 
x°  —  x")  ,  and  d  is  the  spatial  dimension.  The  RK  shape  functions  are 

constructed  using  monomials  as  basis  functions  expressed  as 


4(x)  =  <pQ(x-x7) 


\a\<n 


(85) 


where  <pa(x— x7)  is  the  kernel  function  with  support  size  a,  {(x  — x/)a}|a<n 

is  the  set  of  monomial  basis  functions,  and  {^a(x)}|a|<n  are  the  coefficients 

of  the  basis  functions  obtained  by  imposing  the  following  reproducing 
conditions 


y4(xV=xa,  |a|<n  (86) 

I 

Obtaining  ba(x)  from  Equation  86,  the  RK  shape  function  is  constructed 

4  (x)  =  Hr  (O)ivr1  (x)Hr  (x  -  Xj  )(pa  (x  -  x7 )  (87) 

where  H(x)  =  {xa}a<n  ={l,x1,---,xd}  is  the  vector  containing  all  the  basic 
functions  and  M  is  the  moment  matrix  defined  as 

M(x)  =  y  H(x  -  x,  )Hr  (x  -  Xj  )cpa  (x  -  Xj )  (88) 

I 

To  consider  enrichment  of  RKPM  by  a  set  of  enrichment  functions,  p,  (x) , 
the  following  reproducing  conditions  are  introduced  (Fleming  et  al.  1997) 


Qi  00  =  y  4  OOQi  (x , ) + Pi  (x) 


(89) 
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where  Q,  (x)  is  the  target  function  and  p;(x)  is  the  corresponding 
enrichment.  The  enrichment  is  thus  expressed  as 

Pi(x)  =  Qf(x)-J]4(x)Ql.(xJ)  (90) 

I 

The  enriched  RKPM  is  given  in  the  following  form 

u\x)  =  J2tiWui  (91) 

For  the  microscopic  cell  analysis,  the  following  target  functions  for 
enrichment  of  the  crack  tip  solution  (Moes  and  Belytschko  2002)  are 
considered 


[Q„Q2,QM  = 
[Qi>Q2’Q3]  = 


\fr  sin  ^  ,yfr  cos  ^  ,Vr  sin  ^  sin0,\/rcos 


sin0 


for  T.F.FM 


(92) 


6 

3 

6 

0 

°) 

rsin 

,2, 

,r2  sin 

,2, 

,r  sin 

A 

for  Cohesive  Crack  Model 


The  functional  for  the  microscopic  cell  problem  is  expressed  as 

n  =  |/‘  Cg :Cijklskldn  -  £  utttdr + 1 £  (u,  -  u, )(u,  -  iZ,  )dr  (93) 

where  «  is  the  penalty  parameter  for  imposing  the  essential  boundary 
conditions.  Introducing  the  enriched  RK  approximation  in  Equation  91 
into  the  stationary  condition  of  Equation93,  we  have 

(K  +  aKu)U  =  f  +  afu  (94) 

where  K  is  the  stiffness  matrix,  f  is  the  force  vector,  and  Ku  and  fu  are 

the  terms  associated  with  imposition  of  essential  boundary  conditions.  For 
the  simulation  of  solids,  the  penalty  parameter  a  is  chosen  to  be 
103  — 106  E  where  E  is  the  elastic  modulus. 

Numerical  examples 

In  the  aforementioned  homogenization  calculation  of  damage  parameters, 
modeling  of  the  microscopic  cell  with  microcracks  propagation  is  first 
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performed.  The  microscopic  cell  analysis  results  are  then  processed 
according  to  the  given  homogenization  procedures  to  obtain  the  evolution 
of  damage  parameters.  It  is  noted  that  the  homogenized  damage  evolution 
functions  need  to  be  properly  scaled  with  respect  to  the  ratio  between 
dimensions  of  the  microscopic  cell  and  mesh  to  avoid  the  “numerical  size 
effect.”  This  is  discussed  in  the  following  section. 

A  notched  beam  under  three  point  bending,  as  shown  in  Figure  16,  is  to  be 
modeled  by  the  proposed  methods.  This  problem  was  extensively  studied 
both  experimentally  (Petersson  1981)  and  numerically  (de  Borst  1986; 
Meyer  et  al.  1994).  The  overall  behavior  of  the  beam  is  governed  by  the 
Mode  I  crack  initiated  at  the  vertex  of  the  notch.  Hence  the  microscopic  cell 
model  with  an  evolving  Mode  I  crack  is  considered  for  generation  of  the 
damage  evolution  function.  In  this  example,  we  consider  the  two-parameter 
damage  model  as  previously  given  with  the  compressive  damage  evolution 
suppressed  due  to  the  bending  condition  and  moderate  beam  thickness.  The 
material  properties  for  elastic  modulus  and  Poisson’s  ratio  are  30  GPa  and 
0.2,  respectively.  Tensile  strength  of  the  material  is  fu  =3.33  MPa,  and 

fracture  energy  of  the  cohesive  crack  is  G,  =  124  N/m  (de  Borst  1986). 


Characterization  of  damage  evolution  functions 

The  microscopic  cell  geometry  with  initial  crack  and  boundary  conditions 
are  shown  in  Figure  17.  We  consider  a  microscopic  cell  subjected  to  uniform 
tension  with  a  centered-crack  propagating  perpendicular  to  the  loading 
direction.  In  this  study,  the  initial  crack  length  is  set  to  be  one-tenth  of  the 
microscopic  cell  dimension.  Based  on  the  numerical  method  described 
earlier,  we  first  obtain  the  stress  and  strain  fields  in  the  microscopic  cell  at 
different  loading  stages.  Normal  stresses  in  the  direction  of  loading  at 
different  stages  of  crack  propagation  are  shown  in  Figure  18.  The 
homogenized  stress  and  strain  are  shown  in  Figure  19(a). 
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Figure  17.  Microscopic  cell  problem  and  cohesive  crack  model. 


(a)  (b)  (c)  (d)  (e) 

Figure  18.  Normal  stress  contours  (in  direction  of  loading)  at  different  loading  stages. 


Figure  19.  (a)  Homogenized  stress  and  strain  and  (b)  tensile  damage  parameter 
calculated  by  the  microscopic  stress  and  strain  fields,  where  a,  b,  c,  d  and 
e  denote  the  corresponding  homogenized  stresses,  strains  and  damage 
parameters  calculated  based  on  the  stress  fields  (a),  (b),  (c),  (d)  and 
(e)  shown  in  Figure  18,  respectively. 
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The  energy  bridging  Equation  69  is  employed  to  calculate  the  damage 
variable  as  a  function  of  deformation  as  previously  discussed.  The 
homogenized  tensile  damage  evolution  is  shown  in  Figure  19(b). 

Here  we  study  the  size  effect  by  considering  various  sizes  of  microscopic 
cells  which  correspond  to  structural  models  with  coarse,  medium  and  fine 
meshes.  For  this  purpose,  define  a  dimensionless  parameter,  A ,  as  the 
ratio  between  the  microscopic  length  parameter,  lmic ,  and  the  macroscopic 

length  parameter,  lmc .  Here  we  consider  the  dimension  of  the  microscopic 
cell,  / ,  and  height  of  the  beam,  h ,  as  the  two  scale  parameters  and  define 
their  ratio  as 


A 


^ mic 


l 

h 


(95) 


It  is  observed  in  Figure  20  that  the  homogenized  stress-strain  curves  are 
strongly  affected  by  the  size  of  the  microscopic  cell.  The  cohesive  law 
employed  in  the  microscopic  analysis  of  crack  propagation  involves  a  length 
scale,  which  is  the  crack  opening  displacement  that  corresponds  to  zero 
stress  ( 2Gj  If  in  Figure  17),  called  the  critical  crack  opening  displacement. 

This  length  scale  does  not  scale  with  the  microscopic  cell,  and  thus  leads  to 
different  homogenized  stress-strain  curves  based  on  different  cell  dimen¬ 
sions.  As  the  cell  dimension  increases,  the  homogenized  strain  cor¬ 
responding  to  the  critical  crack  opening  displacement  decreases,  and  yields 
strain  softening  with  a  larger  negative  slope  as  shown  in  Figure  20. 
Computationally,  the  macroscopic  mesh  dimension  is  used  to  represent  the 
averaged  material  behavior  within  the  mesh.  Thus,  if  a  cohesive  law  is  used 
in  a  microscopic  cell  for  obtaining  the  homogenized  stress-strain  curve  for  a 
macroscopic  calculation,  the  corresponding  microscopic  cell  dimension 
needs  to  be  dimensionally  equivalent  to  the  mesh  dimension  when  strain 
softening  exists.  Unfortunately  this  is  practically  tedious  for  arbitrary  mesh 
geometry.  Therefore,  a  scaling  law  will  be  introduced  in  the  following 
section  so  that  the  homogenized  stress-strain  curves  for  different  mesh 
points  are  scaled  based  on  a  “reference  microscopic  cell  analysis.” 
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Figure  20.  Homogenized  stress-strain  curve. 


It  is  also  noted  that  as  the  microscopic  cell  size  increases,  the  energy 
dissipation  capacity  of  the  microscopic  cell  (the  area  under  the  stress-strain 
curve)  decreases.  The  reason  is  that  the  energy  dissipated  by  the  cohesive 
crack  propagation  is  dominated  by  the  elastic  energy  within  the  microscopic 
cell  as  its  size  increases.  For  the  microscopic  cell  with  a  dimension  excee¬ 
ding  the  dimension  of  the  macroscopic  solid,  its  elastic  energy  becomes 
greater  than  the  cohesive  crack  opening  energy,  hence  the  snap-back 
behavior  is  observed  as  shown  in  Figure  20  for  A  =  5,10.  Figure  21(a) 

demonstrates  the  size  dependence  of  the  calculated  nominal  strength  that 
can  be  well  fitted  to  the  size  effect  law  proposed  in  (Bazant  1984) 


°N  = 


B  =  - 

JiTTP  k 


where  fu  is  the  tensile  strength  of  concrete,  /  is  the  specimen  dimension, 
and  B  and  /„  are  material  parameters  identified  by  experimental  data  or 

numerical  simulation.  This  size-dependent  property  is  due  to  influence  of 
the  internal  length  scale,  i.e.,  the  crack  opening  displacement  characteristic 
length,  which  does  not  scale  with  the  overall  dimension  of  the  microscopic 
cell  and  specimen.  The  scale  dependence  of  rupture  strain,  su ,  is  also  shown 

in  Figure  21(b),  where  a  size  effect  law  for  su  can  also  be  extracted  from  the 
numerical  results. 
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Figure  21.  Size  dependence  of  microscopic  cells:  (a)  nominal  strength; 
(b)  difference  between  rupture  strain  eu  and  peak  strain  et. 


For  microscopic  cells  with  different  dimensions  and  with  the  center  crack 
dimension  defined  in  proportion  with  the  microscopic  cell  dimension,  the 
homogenized  damage  evolution  curves  takes  different  paths  from  o  to  1,  as 
can  be  seen  in  Figure  22.  Note  that  a  critical  size  lc  for  the  microscopic  cell 

exists  (Axo.05  in  Figure  22).  This  inherent  critical  size  represents  the 
constraint  on  the  discretization  of  the  macroscopic  structure. 


Figure  22.  Microcracks  informed  damage  evolutions. 
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Remark  3.4 

Besides  the  dimension  of  the  microscopic  cell,  the  characterized 
macroscopic  damage  evolution  is  also  affected  by  the  criteria  of  the 
microcrack  propagation.  While  cohesive  energy  was  employed  for  this 
work  as  the  crack  propagation  criterion,  for  general  applications  it  should 
be  carefully  investigated  according  to  the  key  characteristics  of  material 
behavior. 

Mesh  insensitive  solution  by  scaled  damage  evolution  functions 

The  characterized  tensile  damage  evolution  equation  for  the  two-parameter 
damage  model  is  employed  in  the  structural  level  analysis  based  on 
continuum  damage  mechanics.  Due  to  softening  behavior  of  the  material 
response,  the  arc-length  method  is  used  for  the  nonlinear  analysis.  To  study 
mesh  sensitivity  of  the  multiscale  analysis,  coarse,  medium  and  fine  meshes 
are  constructed  for  the  notched  beam  as  shown  in  Figure  23. 


Figure  23.  Finite  element  mesh  with  different  levels  of  refinement:  (a)  coarse  mesh; 
(b)  medium  mesh;  (c)  fine  mesh. 


In  the  conventional  damage  models,  the  damage  evolution  curves  are 
directly  used  in  the  structural  analysis  without  consideration  of  the  relation¬ 
ship  between  microstructure  dimensions  and  mesh  size.  Figure  24  shows  a 
strong  mesh  dependency  induced  by  this  standard  procedure  where  only 
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Mid-span  displacement  (m) 


Figure  24.  Mesh  dependent  load-displacement  responses  using  inconsistent  microscopic  cell. 

one  microscopic  cell  is  used  to  characterize  the  damage  evolution  equation 
without  use  of  the  scaling  law  in  Equation  96  and  linear  fit  of  Figure  21(b). 
With  the  proposed  method,  the  microcracks  induced  damage  evolution 
curve  is  first  characterized  by  the  cracked  microscopic  cell  simulation 
results.  By  introducing  the  scaled  damage  evolution  curves  in  Figure  22 
according  to  the  mesh  dimension,  the  mesh  independent  results  are 
obtained.  The  agreement  between  numerical  results  with  different  mesh 
refinements  and  the  benchmark  results  (Abaqus  2005)  is  shown  in 
Figure  25.  The  stress  and  strain  contours  representing  damage  evolution  in 
the  structures  are  shown  in  Figure  26  and  Figure  27. 

Remark  3.5 

The  macroscopic  mesh  is  a  numerical  representation  of  the  solids  included 
within  its  domain,  and  the  stress-strain  relation  introduced  in  the 
quadrature  points  describes  the  homogenized  behavior  of  the  solid  within 
the  mesh.  For  an  elastic  solid,  the  homogenized  stress-strain  curve  is 
insensitive  to  the  mesh  dimension,  and  thus,  no  mesh  dependency  issues 
exist  for  elastic  problems.  On  the  other  hand,  the  softening  solid  typically 
involves  length  scales  in  the  material  laws,  for  example,  the  ones 
constructed  by  homogenization  of  microscopic  fracture  analysis  with 
cohesive  law  introduced  in  this  work.  When  the  length  scales  in  the  material 
laws  do  not  scale  properly  with  the  representative  domain  of  the 
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macroscopic  quadrature  points,  mesh  dependent  results  arise.  In  the 
present  work,  the  mesh  dependency  in  the  softening  problem  is  removed  by 
adopting  a  scaling  law  to  a  “reference  homogenized  stress-strain  curve” 
obtained  from  a  reference  microscopic  cell  analysis. 


Mid-span  displacement  (m) 

Figure  25.  Mesh  independent  load-displacement  responses  using  consistent 

microscopic  cell. 
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Figure  26.  Contour  of  normal  stress  in  horizontal  direction. 


LE,  LE11 

(Ave.  Crit.:  75%) 
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Figure  27.  Contour  of  normal  strain  in  horizontal  direction. 
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4  Enhanced  Frictional  Natural  Kernel 
Contact  Algorithm  for  Impact  Modeling 

Kernel  contact  algorithm 

Conventional  contact  algorithms,  such  as  the  penalty  method,  require  the 
potential  contact  surfaces  to  be  pre-defined.  However,  for  penetration 
problems,  contact  surfaces  are  changing  continuously  and  they  cannot 
always  be  defined  a  priori.  In  this  study,  a  kernel  contact  algorithm  was 
introduced  that  utilizes  the  interaction  of  kernel  functions  between 
contacting  bodies  to  naturally  serve  as  the  impenetration  condition. 

Consider  a  semi-Lagrangian  discretization  of  a  continuum  by  a  set  of 
RKPM  points,  as  shown  in  Figure  28,  with  each  point  carrying  nodal 
volume  Aj ,  mass  m7 ,  kernel  function  cpa  (x  -  x1 ) ,  and  state  and  field 

variables.  The  interaction  between  RKPM  points  (Figure  28c)  via  the 
overlap  of  kernel  supports  induces  a  stress 

d(x)  =  ^DB/(x)d/  (97) 

I 

where  D  is  an  artificial  material  response  tensor  introduced  between  two 
contacting  bodies,  and  B7  (x)  is  the  smoothed  gradient  of  the  shape  function. 

The  layer  of  the  artificial  material  serves  as  a  medium  to  transmit  the 
contact  forces. 

An  elastic-perfect  plastic  material  model  in  the  contact  layer  is  introduced 
to  represent  the  frictional  contact  condition  without  imposing  kinematic 
constraints.  The  yield  function  is  defined  as 

/(t)  =  M-/i|611|  (98) 

where  t  =  [d  12  a  13  ] ,  ||x||  =  (x  •  x)±/2 ,  and  a  is  the  projected  Cauchy  stress 

tensor  onto  the  local  coordinate  in  which  the  1 -direction  is  in  alignment 
with  the  outward  normal  of  the  contact  surface,  and  2-  and  3-  directions 
are  in  alignment  with  two  mutually  orthogonal  vectors  lying  in  the  contact 
surface  (Figure  28d). 
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fa)  Two  sapordtad  bodies  in  motion 


<(b}  Kernel  overlapping 


fd)  Surface  normal  and  tangential 
direclions  in  contact  processing  zoni 


(t)  Keffiel  cortaot  as  impanelalion 
condition 


Figure  28.  Kernel  contact  algorithm  by  kernel  interaction  between  contacting  bodies. 


As  such,  the  yield  stress  as  /n\an\  mimics  the  friction  stress  induced  by 
normal  stress  <rn ,  and  the  sliding  condition  is  represented  by  the  yield 
condition  in  the  plastic  model. 


stick  condition  (elastic) 

t|| I  sliding  condition  (plastic) 


With  isotropic  hardening  and  perfect  plasticity  assumptions,  this  plasticity 
model  represents  the  Coulomb’s  friction  law. 


If  SCNI  or  SNNI  is  introduced  for  the  domain  integration,  the  internal 
force  acting  on  point  I  is 


f/=  EB/(xXxjK 


(100) 
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where  JV7  =  |j  |  cpa  (x7  —  x7)  ^  O}  is  the  set  containing  neighbors  of  point  I. 

In  the  previous  equation,  the  total  force  acting  on  point  I  is  obtained  by 
summing  up  all  pair  interactions  between  point  I  and  its  neighbors.  This 
property  is  applied  to  the  interaction  between  contacting  bodies  as  follows 

fj  =  E^(xjMxj)A/>  VJgG7  and  J^G7  if  n(xJ)-d(xJ)-n(xJ)<0  (101) 

JeNj 


where  G7  is  the  group  of  points  discretizing  the  body  that  particle  I 
belongs  to,  and  n(xJ)  is  the  unit  contact  surface  normal.  The  unit  normal 

vector  will  be  defined  by  a  level-set  method  in  the  next  section.  In  this 
approach,  when  two  bodies  are  approaching  each  other,  the  pair-wise 
interactions  due  to  overlapping  kernel  functions  serve  as  a  natural 
impenetration  condition  as  shown  in  Figure  28c.  The  radius  of  kernel 
support  determines  the  numerical  length  scale  in  the  normal  contact.  The 
stick  and  slip  conditions  can  be  calculated  based  on  the  tangential  stress 
n  at  in  the  contact  processing  zone.  Since  the  contacting  bodies  exhibit 
high  velocity  gradients  across  the  contact  surface,  stability  conditions 
given  in  Chapter  2  are  crucial  for  the  kernel  contact  algorithm  to  be  stable. 

Level  set  algorithm  for  determination  of  surface  normal 

The  determination  of  contact  surface  using  point-based  discretization  of 
contacting  bodies  in  the  RKPM  requires  additional  effort.  The  accuracy  of 
surface  normal  estimation  at  discrete  points  based  on  nodal  data  is  crucial 
for  contact  force  calculation.  In  this  section,  a  level  set  method  is  introduced 
to  estimate  the  contact  surface  normal  under  the  RK-based  kernel  contact 
framework. 

Consider  the  discretization  of  a  level  set  function 

^(X)  =  E^(X)C/’  f  £  G7  U G2  (102) 

I 

where  GK  is  the  group  of  points  contained  in  body  K,  and  C7  is  the  level 
set  nodal  value  associated  with  the  RK  shape  function  'F/  (Figure  29a). 
The  level  set  nodal  value  is  defined  as 
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^  :  points  in  G* 

(§)  (•)  :  contact  points  - 

-►  :  outward  normal 

^  :  points  in  G2 

:  contact  region 

Figure  29.  Level  set  algorithm  to  identify  contact  nodes  and  obtain  normal  vector  of  the 

contact  surface. 


[  1  ifIeG1 
{-1  ifIeG2 


(103) 


The  level  set  function  in  Equation  102  gives  a  zero  level  set  between  these 
two  bodies,  which  severs  as  the  contact  surface,  as  can  be  seen  in  Figure  29b 
and  29c.  The  contact  surface  outward  normal,  shown  in  Figure  29c,  then 
can  be  estimated  by 


n  =  -V^(x)/||v^(x)||  (104) 

Note  that  the  gradient  operator  in  Equation  104  can  be  replaced  by  the 
smoothed  gradient  operator  described  in  Chapter  2.  The  contact  force  can 
be  therefore  applied  to  the  contact  points  following  the  frictional  kernel 
contact  algorithm  described  in  the  previous  section. 

Numerical  examples 

Several  numerical  examples  solved  by  the  frictional  natural  kernel  contact 
algorithm  are  presented  in  this  section  to  demonstrate  the  accuracy  of  the 
proposed  method. 
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Elastic  collision  of  two  bodies 

Consider  a  two-body  elastic  contact  problem  as  illustrated  in  Figure  30.  The 
two  bodies  have  the  same  discretizations  and  same  material  properties 
(Young’s  modulus  and  mass  density).  Body  1  is  prescribed  with  initial 
velocity  of  500  m/sec  and  Body  2  is  initially  at  rest.  First  of  all,  the  natural 
kernel  contact  algorithm  is  examined  by  the  case  with  same  support  sizes  of 
the  two  bodies.  Numerical  results  of  the  two-body  elastic  collision  are  pre¬ 
sented  in  Figure  31.  As  observed  from  the  results,  the  contact  of  two  bodies 
is  handled  properly  with  the  proposed  natural  kernel  contact  algorithm, 
whereas  both  linear  momentum  and  total  energy  are  conserved  during  the 
contact  process. 


Figure  31.  Numerical  results  of  two-body  elastic  contact  with  same 
properties  and  support  sizes. 
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To  demonstrate  the  effect  of  tensile  contact  force  correction  which  is  the 
condition  as  indicated  in  Equation  101,  another  two-body  contact  problem 
with  masses  m^m  and  m2  =  4m ,  and  normalized  support  sizes  ax  =  2.0 

and  a2  =1.75  is  investigated.  Figure  32  shows  the  numerical  results  with 

and  without  correction  of  tensile  contact  force.  It  is  clear  that  the  pulling 
back  effect  due  to  non-physical  tensile  contact  force  between  two  bodies 
can  be  eliminated  by  introducing  the  tensile  force  correction,  and  the  post¬ 
contact  velocities  of  two  bodies  are  close  to  the  theoretical  values. 


(a)  Without  correction  of  tensile  contact  force  (b)  With  correction  of  tensile  contact  force 
V,  =  -265.2  m/sec,  Vj  =  1 91.3  ra/sec  Vt  =  -299.9  ni/sec,  V3  =  192.6  m/sec 


Figure  32.  Numerical  results  of  two-body  elastic  contact  with 

mx  =m,  m2  =4 m  and  ax  =2.0,  a2  =1.75 

Sliding  block  on  a  rigid  surface 

An  elastic  block  resting  on  a  rigid  surface  is  illustrated  in  Figure  33.  Both 
the  rigid  surface  and  flexible  body  are  discretized  by  RKPM  particles.  After 
the  block  rests  on  the  rigid  wall  and  reaches  equilibrium,  the  rigid  surface 
rotates  an  angle  of  60  deg  and  the  block  begins  to  slide.  Frictional 
coefficients  jU  =  0.0  and  jU  =  0.2  are  considered.  Figure  34  shows  the 
numerical  results  obtained  by  RKPM  simulation  for  the  case  ju  =  0.2 ,  and 

the  numerical  predictions  of  the  displacement  of  the  elastic  body  compared 
with  analytical  solutions  are  illustrated  in  Figure  35  for  both  //  =  0.0  and 
//  =  0.2 .  It  is  observed  that  the  numerical  solutions  have  good  agreement 
with  the  analytical  solutions. 
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Figure  33.  Schematic  of  an  elastic  block  resting  on  a  rigid  surface 
subject  to  gravity. 


Figure  34.  Sliding  block  problem  with  the  enhanced  natural  kernel  contact  algorithm  and  level 
set  identification  of  contact  surface  normal  ( //  =  0.2 ). 
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Figure  35.  Comparison  of  block  displacement  in  the  sliding  block  problem  with  the  natural 
kernel  contact  algorithm  and  level  set  identification  of  contact  surface  normal. 


Dynamic  Hertz  contact 

The  dynamic  Hertz  contact  problem  is  analyzed  with  an  elastic  ball 
impacting  toward  a  rigid  surface  as  illustrated  in  Figure  36.  Figure  36  also 
shows  the  RKPM  discretization  of  the  elastic  ball  and  rigid  wall.  The 
numerical  solution  of  contact  radius  is  compared  with  analytical  solution 
(Timoshenko  and  Goodier  1934)  as  shown  in  Figure  37.  As  seen  from  the 
results,  the  proposed  kernel  contact  algorithm  with  contact  surface 
identification  via  level  set  method  is  effective  for  the  dynamic  contacts. 


Figure  36.  Hertz  contact  problem:  (a)  problem  statement  and  (b)  RKPM  discretization. 
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Figure  37.  Numerical  result  of  contact  radius  compared  to 
analytical  solution  (Timoshenko  and  Goodier  1934). 


Dynamic  contact  of  elastic  ring 

Analyzed  here  is  an  elastic  ring  with  Young’s  modulus  E  =  10  GPa  and 
Poisson’s  ratio  v  =  0.4  impacting  on  a  rigid  surface  with  the  initial  velocity 
v0  =  50  m/sec  as  depicted  in  Figure  38.  Both  gravitational  forces  and 

frictions  are  ignored  in  the  simulation.  The  progressive  impact  processes 
for  the  contact  surface  normal  determined  by  level  set  algorithm  are 
presented  in  Figure  39.  Figure  40  shows  the  comparison  of  the  impact  and 
reflected  angles  between  the  two  algorithms  of  surface  normal  calculation, 
one  by  nodal  orientation  between  two  particles  and  the  other  one  by  level 
set  algorithm  as  introduced  in  the  previous  section.  It  is  clear  that  the  sim¬ 
ple  computation  of  contact  surface  normal  by  level  set  algorithm  provides 
fairly  accurate  surface  normal  for  impact  modeling. 

Next,  the  dynamic  contact  of  two  elastic  rings  is  analyzed  as  shown  in 
Figure  41.  The  rings  are  made  of  the  same  material  properties  as  the 
previous  one-ring  problem.  The  two  rings  are  moving  toward  each  other 
with  an  initial  velocity  v0  =  50  m/ sec .  The  progressive  contact  processes 

and  total  energy  of  the  system  are  illustrated  in  Figures  42  and  43, 
respectively.  It  is  seen  that  the  contact  of  these  two  elastic  rings  is  properly 
handled  with  the  proposed  kernel  contact  algorithm. 
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Figure  39.  Impact  of  an  elastic  ring  on  a  rigid  surface  by  natural  kernel  contact 

algorithm. 
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(a)  Surface  norm  a  l  obtained  by  nodal 
orientation  between  two  particles 


(l>|  Surface  normal  computed 
by  level  set  algorithm 


Figure  40.  Comparison  of  impact  and  reflected  angles  between  two  algorithms  for 
determination  of  contact  surface  normal. 


Figure  42.  Impact  of  two  elastic  rings  using  natural  kernel  contact  algorithm. 
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Figure  43.  Time-history  of  total  energy  for  the  two-ring  impact  problem. 

Cylindrical  aluminum  bar  impacting  a  rigid  wall 

This  classical  impact  problem  with  available  experimental  and  numerical 
results  (Taylor  1948;  Wilkins  and  Guinan  1973)  is  used  to  test  the 
performance  of  the  proposed  contact  algorithm  for  plastic  deformation. 
The  initial  radius  and  initial  height  of  the  aluminum  cylindrical  bar  are 
0.391  and  2.346  cm,  respectively.  The  material  properties  of  the  cylinder 
are:  density  p  =  2700  kg/m3 ,  Young’s  modulus  E  —  78.2  GPa ,  Poisson’s 
ratio  v  =  0.3  ,  and  J2  plasticity  with  yield  strength  oY  =  0.29  GPa .  The 

initial  impact  velocity  is  373  m/sec  and  the  rigid  surface  is  assumed  to  be 
frictionless.  Perfect  plasticity  and  strain-hardening  elasto-plasticity  with 
the  following  hardening  rules  are  considered  for  this  problem 

H(ep)  =  0  (105) 

K(ep)  =  aT(l  +  125ep)°'01  (106) 


where  ep  is  the  effective  plastic  strain,  and  H  and  K  are  the  plastic 
modulus  and  yield  stress,  respectively  (Chen  et  al.  1996). 
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Both  the  aluminum  bar  and  rigid  plate  are  discretized  by  the  RKPM 
particles  in  three-dimensions  (29,788  RKPM  particles)  as  illustrated  in 
Figure  44.  The  contact  between  the  aluminum  bar  and  the  rigid  wall  is 
handled  by  the  frictional  natural  kernel  contact  algorithm,  and  the 
problem  is  modeled  by  the  semi-Lagrangian  formulation  with  SNNI  nodal 
integration. 


Figure  44.  Schematic  of  Taylor  bar  problem  and  corresponding  RKPM  discretization. 


Table  1  summarizes  the  comparison  of  deformed  heights  and  radii 
between  RKPM  predictions  and  experimental  measurement.  The 
deformed  configurations  of  the  aluminum  bar  at  different  time  frames  for 
both  hardening  and  perfect  plasticity  materials  obtained  by  RKPM  SNNI 
semi-Lagrangian  calculation  are  shown  in  Figure  45.  This  study  also 
compares  Lagrangian  RKPM  results  reported  by  Chen  et  al.  (1996).  Note 
that  for  problems  with  moderate  level  of  deformation,  such  as  in  this 
impact  problem,  the  semi-Lagrangian  method  is  at  best  as  good  as  the 
Lagrangian  method.  The  results  in  Table  1  show  that  the  semi-Lagrangian 
RKPM  results  are  very  close  to  the  Lagrangian  RKPM  results,  indicating 
the  robustness  of  the  proposed  semi-Lagrangian  approach. 


ERDC/GSL  TR-11-35 


53 


Table  1.  Comparison  of  deformed  geometries  for  Taylor  bar  impact  problem. 


Perfect  Plasticity 

Hardening 

RKPM  Lagrangian 
(Chen  1996)* 

RKPM  Semi- 
Lagrangian 

SNNI 

RKPM  Lagrangian 
(Chen  1996) a 

RKPM  Semi- 
Lagrangian 

SNNI 

Experiment 
(Wilkins  and 

Guinan  1973) 

Height  (cm) 

1.454 

1.442 

1.645 

1.642 

1.651 

Radius  (cm) 

1.051 

1.025 

0.837 

0.819 

NA 

a  Axisymmetric  model  with  11x31  particles  was  analyzed  in  (Chen  1996). 


(a)  Hardening 
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X  X  X  X 


(b)  Perfect  plasticity 


Figure  45.  Cylindrical  impact  deformations  predicted  by  the  semi-Lagrangian  SNNI  RKPM. 
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5  V&V  Penetration  Study 

An  extensive  verification  and  validation  (V&V)  study  was  conducted  to 
validate  performance  of  the  Nonlinear  Meshfree  Analysis  Program 
(NMAP)  parallel  code  (Chi  et  al.,  in  preparation)  for  penetration  modeling. 
The  V&V  problem  set  consisted  of  18  simulations  that  followed 
penetration  experiments  (Boone,  in  preparation)  conducted  in  ERDC’s 
Fragment  and  Small  Arms  laboratory.  Results  from  two  experiments  were 
provided  by  ERDC  at  the  beginning  of  the  V&V  study,  and  results  from  five 
additional  tests  were  provided  approximately  half-way  through  the  V&V 
effort.  The  remaining  n  experiments  were  modeled  blindly.  Calculations 
were  jointly  performed  on  the  UCLA  Hoffman2  cluster  and  systems  at  the 
ERDC  DoD  Super  Computing  Resource  Center. 

Experimental  setup 

The  penetration  experiments  were  conducted  by  firing  a  spherical 
projectile  into  a  thin-panel  concrete  target  with  an  impact  angle  of  zero- 
degree  obliquity.  Impact  velocity,  projectile  size,  and  panel  thickness  were 
varied  to  obtain  terminal  ballistic  conditions  ranging  from  strong 
overmatch  (i.e.,  projectile  exited  the  panel  with  significant  residual 
velocity)  to  stopping  the  projectile.  Three  steel  projectile  sizes  were  used, 
with  diameters  of  7.9  mm  (32  grain  mass),  11.1  mm  (86  grain  mass),  and 
12.7  mm  (129  grain  mass).  All  projectiles  were  made  from  American 
Society  for  Testing  and  Materials  (ASTM)  A681  (2008),  Grade  200 
hardened  S2  tool  steel  with  material  properties  provided  by  ERDC  and 
given  in  Table  2. 


Table  2.  Spherical  steel  projectile  properties. 


Property 

Value3 

Tensile  strength,  ultimate 

2,150  MPa 

Tensile  strength,  yield 

2,000  MPa 

Elongation  at  break 

7%  in  50  mm 

Bulk  modulus 

140  GPa 

Shear  modulus 

80  GPa 

a  Material  property  data  provided  by  ERDC. 


The  target  panels  were  nominally  305  mm  x  305  mm  in  planimetric 
dimension,  with  thicknesses  of  approximately  12.7  mm,  25.4  mm  or 
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38.1  mm.  All  panels  were  produced  using  CorTuf  ultra  high-strength 
concrete  (Williams  et  al.  2009).  The  CorTuf  material  was  reinforced  with 
approximately  30-mm  long,  hooked-end  Bekaert  Dramix®  ZP305  steel 
fibers.  Due  to  the  length  of  the  steel  fibers  with  respect  to  the  panel 
thicknesses,  a  certain  amount  of  preferential  fiber  orientation  (in  the  plane 
of  the  panels)  occurred.  However,  the  degree  of  preferential  orientation 
was  not  quantified. 

In  each  experiment  the  projectile  was  fired  from  an  un-rifled  barrel; 
therefore  no  rotational  velocity  was  induced.  The  target  panels  were 
mounted  in  a  support  fixture  and  were  held  in  place  by  clamping  all  four 
edges  between  a  pair  of  metal  plates.  Rubber  strips  were  placed  between  the 
plates  and  panel  to  avoid  stress  concentrations  resulting  from  the  uneven 
concrete  surfaces.  The  grip  of  the  plates  onto  the  panels  was  approximately 
25  mm  on  all  sides.  With  the  exception  of  Trial  2  and  Run  15,  projectiles 
were  fired  at  the  center  of  the  CurTuf  targets.  In  Trial  2  and  Run  15,  the 
projectiles  were  fired  to  impact  the  panels  in  the  top  right  and  bottom  right 
quadrants,  respectively.  In  both  cases  the  projectile  impacted  approximately 
75  mm  from  the  panel  edges.  Quantification  of  the  experimental  results 
included  measurement  of  the  impact  and  exit  velocities,  impact  and  exit 
crater  sizes,  tunnel  size,  and  post-test  projectile  mass. 

A  CorTuf  target  panel  mounted  into  the  test  support  fixture  is  shown  in 
Figure  46. 


CorTuf  panel 


metal  plate 


X 


(a)  (b) 

Figure  46.  CorTuf  panel  mounted  in  test  support  fixture,  (a)  front  view 
(impact  face)  and  (b)  side  view. 
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NMAP  model 

To  construct  the  NMAP  models,  the  12.7-mm-,  25.4-mm-,  and  38.1-mm- 
thick  panels  were  discretized  using  107,811;  190,000;  and  357,309  nodes, 
respectively.  Nodal  spacing  in  the  plan  dimensions  was  increased  from  the 
center  toward  the  outer  edges,  with  a  maximum  change  in  nodal  spacing  of 
4%  between  each  node.  Nodal  spacing  was  held  constant  through  the  panel 
thickness  as  summarized  in  Table  8.  The  projectiles  were  discretized  using 
561;  1,163;  and  1,517  nodes  for  the  32-grain,  86-grain,  and  129-grain 
projectiles,  respectively.  The  optimum  discretization  for  each  panel  and 
projectile  size  was  studied,  which  is  discussed  in  the  following  section. 
Boundary  conditions  were  applied  to  the  panels  by  fixing  nodal  displace¬ 
ments  along  all  panel  edges. 

The  projectiles  were  modeled  using  a  J2  plasticity  model  with  isotropic 
hardening.  Material  yield  stress  was  modified  to  account  for  the  so-called 
dynamic  increase  factor  (DIF)  associated  with  strain-rate  effects.  A  DIF  of 
either  1.2  or  1.5  was  used  depending  on  projectile  impact  velocity,  which  is 
discussed  in  the  following  section.  Material  model  parameters  used  for  the 
projectile  are  given  in  Table  3.  The  CorTuf  panels  were  modeled  using  the 
MIDM  (microcrack  informed  damage  model)-enhanced  AFC  model  (Ren 
et  al.  2011),  with  selected  experiments  also  modeled  with  the  original  AFC 
model  (Adley  et  al.  2010).  AFC  model  parameters  for  the  CorTuf  material 
were  provided  by  ERDC  and  are  given  in  Table  4. 


Table  3.  J2  material  model  parameters  for  projectile. 


Parameter 

Value 

Young’s  modulus,  E 

200  GPa 

Poisson’s  ratio,  v 

0.26 

Yield  stress,  aya 

2400  MPa/3000  MPa 

Hardening  modulus,  H 

2500  MPa 

Density 

7806  kg/m3 

a  DIF=1.2/DIF=1.5. 


Determination  of  optimal  parameters  for  V&V  modeling 

The  following  parametric  study  was  conducted  to  determine  optimum 
parameters  for  modeling  the  V&V  penetration  problems.  This  study  was 
carried  out  in  order  to  provide  guidelines  for  the  user  to  efficiently 
perform  numerical  analysis  of  impact  problems.  Numerical  studies  are 
provided  to  illustrate  the  choice  of  optimal  parameters. 
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Table  4.  AFC  model  parameters  for  CorTuf  panel. 


Parameter 

Value 

G  (shear  modulus) 

18.457  GPa 

Cl 

1016.3  MPa 

C2 

908.65  MPa 

C3 

0.0125 

C4 

0.10382 

C5 

792.89  MPa 

Q1 

artificial  bulk  viscosity  not  used 

Q2 

artificial  bulk  viscosity  not  used 

PMIN 

6.8947  MPa 

C6 

172.37  MPa 

C7 

0.00781 

C 

7919.2  MPa 

D 

-29.205  GPa 

s 

187.10  GPa 

C9 

77.958  GPa 

CIO 

0.24863 

D1 

4.0597x10-1°  Pa1 

AN 

1.7345x10-9  Pa1 

TXETXCR 

0.625 

PRECRIT 

0.177  xlO22 

HMIN 

1 

RHO  (density)3 

2267.4  kg/m3 

DAMP 

0.0001 

FC  (tensile  strength  for  MIDM) 

6.8947  MPa 

a  2267.4  kg/m3  was  value  used  in  model;  density  reported  by  ERDC  was  2557  kg/m3. 


•  Discretization  of  the  two  bodies: 


Discretization  and  support  size  discrepancies  between  the  contacting 
bodies  were  studied  to  determine  the  effect  on  the  numerical  solution.  Two 
studies  are  conducted  on  two  elastic  bodies  impacting  each  other.  Body  l 
of  mass  rrh  has  initial  velocity  vt.  Body  2  of  mass  m2  is  initially  at  rest. 
Assuming  an  initial  velocity,  Vi,  of  500  m/s,  analytical  results  for  the 
velocity  of  the  bodies  after  contact  can  be  obtained  by 


_  (m1  -m2)v1  +  2m2v2 
1  m1  +  m2 


-300 


(107) 
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_  (m1  -m2)v2  +2m2v1 
2  m1+m2 


200 


(108) 


The  first  study  investigates  the  impact  of  two  elastic  bodies  with  the  same 
discretization,  as  described  in  Figure  47.  Velocity  history  of  the  two  bodies 
is  shown  in  Figure  48.  This  study  reveals  that  if  the  two  bodies  are 
discretized  by  the  same  nodal  density  and  support  size,  the  numerical 
results  from  the  contact  algorithm  are  consistent  with  the  analytical 
solution.  Nevertheless,  for  a  small  penetrator,  using  the  same  nodal 
density  for  projectile  and  panel  leads  to  a  very  fine  discretization  of  the 
panel  and  is  not  efficient.  Therefore,  a  compromise  is  made  by  selecting  a 
coarser  nodal  density  for  the  panel. 


Body  1  j 
m,  =  m 

v !  =  500 

Body  2 

m7  -  4m 

Figure  47.  Study  1  problem  statement. 


Figure  48.  Velocity  histories  of  Study  1. 
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In  the  second  study,  the  case  of  two  elastic  bodies  with  different 
discretizations  is  modeled  as  illustrated  in  Figure  49.  The  effect  of  support 
size  on  the  numerical  results  is  summarized  in  Table  5.  The  results  from 
this  second  study  reveal  that  similar  support  size  for  both  bodies  should  be 
used  even  when  using  different  discretizations.  Conservation  of  linear 
momentum  for  each  support  size  ratio  is  observed  (see  Figure  50),  but 
conservation  of  energy  is  not  satisfied  for  large  support  discrepancies  as 
shown  in  Figure  51. 


Table  5.  Effect  of  support  size  on  numerical  result  in  terms  of  velocity. 


Support  Size3 

Vi 

V2 

2.0  :  1.0 

-152 

163 

2.0  :  1.5 

-294 

199 

2.0  :  1.8 

-299 

200 

a  Support  size  Body  1:  support  size  Body  2. 


Figure  50.  Linear  momentum  history. 
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(a)  support  size  2:1  (b)  support  size  2:1.5  (c)  support  size  2:1.8 

Figure  51.  Energy  history. 


Next,  a  study  was  conducted  with  four  V&V  problems  to  select  the  optimal 
nodal  distance  ratio  between  the  projectile  and  panel.  A  minimum  nodal 
distance  for  the  panel  particles  is  defined  using  the  time-step  obtained 
from  the  Courant-Friedrichs-Lewy  (CFL)  condition  and  the  elastic  and 
plastic  wavelength.  The  one  obtained  considering  the  elastic  wavelength  is 
the  most  conservative.  In  this  discretization,  a  nodal  distance  of 
approximately  l  mm  was  used  for  projectile,  which  gives  a  ratio  of  about 
0.0268  between  nodal  distance  and  the  elastic  wavelength  for  the  system. 
The  numerical  study  in  Table  6  shows  that  the  suitable  projectile  and 
panel  discretization  near  the  contacting  surfaces  corresponds  to  a  nodal 
spacing  in  the  panel  of  approximately  1.4  to  1.5  times  the  nodal  spacing  in 
the  projectile.  The  selected  discretization  of  the  different  projectile  sizes 
according  to  this  ratio  is  shown  in  Table  7. 


Table  6.  Effect  of  discretization  on  exit  velocity. 


Nodal 

Distance  Ratio 

1:1.8a 

1:1.4a 

1:1.2a 

Experiment 

Results 

Number  of 
Particles 

Exit 

Velocity 

(m/sec) 

Number  of 
Particles 

Exit  Velocity 
(m/sec) 

Number  of 
Particles 

Exit  Velocity 
(m/sec) 

Exit  Velocity 
(m/sec) 

Trial  1 

Impact  velocity  = 
1708  m/sec 

183,346 

511.50 

358,472 

571.10 

620,278 

448.00 

684.58 

Trial  2 

Impact  velocity  = 
577  m/sec 

101,019 

12.76 

201,163 

-2.19 

346,763 

5.83 

0 

Run  4 

Impact  velocity  = 
1718  m/sec 

130,892 

570.50 

268,781 

646.00 

519,917 

560.13 

791.56 

Run  5 

Impact  velocity  = 
239.27  m/sec 

59,610 

12.11 

124,157 

3.72 

190,894 

9.78 

0 

a  Projectile  nodal  spacing:panel  nodal  spacing  in  contact  area. 


ERDC/GSL  TR-11-35 


61 


Table  7.  Discretization  of  projectile. 


Projectile  Size 

Nodal  Distance  (m) 

Discretization 

32-grain  -  7.9  mm 

0.0008897 

561  particles 

86-grain  -  11.1  mm 

0.0009698 

1163  particles 

129-grain  -  12.7  mm 

0.0009969 

1517  particles 

•  Time-step  study 

The  stability  analysis  showed  that  the  critical  time-step  size  in  semi- 
Lagrangian  RKPM  is  related  to  the  velocity  gradient,  which  gives  a  time- 
step  approximately  one-tenth  of  that  from  the  CFL  condition  (see  Figure  9). 
This  yields  a  much  more  restrictive  stability  condition  for  contact-impact 
problems  compared  to  the  CFL  condition.  Since  the  stability  condition 
provides  an  estimate  of  the  critical  time-step  based  on  the  result  of  l-D 
VonNeumann  analysis  with  uniform  discretization,  additional  numerical 
experiments  were  performed  on  the  critical  time-step  for  penetration 
modeling.  The  Trial  1  penetration  test  is  used  to  investigate  the  effect  of 
time-step  on  velocity  history,  as  shown  in  Figure  52.  For  the  discretization 
used  in  the  V&V  problems,  the  CFL  time-step  restriction  is  approximately 
2xi0‘7  sec.  Therefore  a  time-step  of  2xio-8  sec  (one  tenth  of  CFL  time-step) 
is  adopted  for  the  low-impact  velocity  cases.  For  the  high-impact  velocity 
cases,  a  smaller  time-step  of  lxicr8  sec  (one-twentieth  of  CFL  time-step)  is 
used  to  capture  the  impact. 

•  Dynamic  Increase  Factor  (DIF) 


Material  yield  stress  of  the  projectile  is  increased  by  a  DIF  to  account  for 
strain-rate  effects.  With  a  lack  of  strain  rate  information  and  experimental 
data  about  the  projectile  material  in  the  V&V  problems,  the  DIF  is  defined 
following  Rule  and  Jones’s  study  (Rule  and  Jones  1998).  This  study 
describes  a  revised  Johnson-Cook  model  and  provides  experimental  data 
for  the  material  behavior  of  steel  at  high  strain  rate,  as  shown  in  Figure  53. 
A  DIF  of  1.2  is  applied  to  the  yield  stress  of  low-impact  velocity  projectiles 
with  velocity  less  than  610  m/sec  and  a  DIF  of  1.5  is  applied  to  the  yield 
stress  of  high-impact  velocity  projectiles  with  velocity  greater  than  610  m/ 
sec.  This  adjustment  of  the  yield  stress  of  the  projectile  is  included  in  input 
file  NMAP_input.DAT.  A  rate-dependent  constitutive  model,  such  as  the 
Johnson-Cook  model,  will  be  considered  in  future  study. 
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Figure  52.  Time-step  effect  on  the  velocity  history. 


Figure  53.  Yield  stress  versus  strain  rate  for  steel  (Rule  and  Jones  1998). 
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•  Boundary  layer  correction: 


The  numerical  contact  surface  implemented  in  NMAP  should  agree  with 
the  physical  contact  surface  in  order  to  define  the  contact  event  as  accurate 
as  possible.  A  correction  of  the  discretization  of  the  two  bodies  is  done  by 
shrinking  half  of  the  length  of  the  boundary  element,  as  illustrated  in 
Figure  54.  The  outer  nodes  of  the  two  bodies  are  moved  inward  with 
adjustment  of  nodal  volume. 


Figure  54.  Boundary  layer  correction. 
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•  Support  size  increase  of  bullet’s  inner  particles: 

Preliminary  studies  of  penetration  simulation  by  NMAP  revealed  a 
nonphysical  separation  of  groups  of  particles  during  projectile  deformation 
as  the  result  of  improper  support  size,  as  shown  in  Figure  55.  This 
phenomenon  is  called  the  kernel  support  instability,  which  is  a  numerical 
instability  due  to  insufficient  kernel  support  coverage  of  neighboring 
particles.  Consider  the  group  of  the  first  layer  of  particles  on  the  surface  of 
the  projectile  to  be  the  “outer  particles,”  and  the  remaining  particles  to  be 
the  group  of  “inner  particles.”  A  remedy  to  the  kernel  support  instability  is 
undertaken  by  increasing  the  kernel  support  size  of  the  inner  particles 
(Figure  56)  in  the  projectile  such  that  it  covers  the  outer  particles.  In  this 
manner,  it  recovers  kernel  support  stability  in  the  projectile  under  high- 
velocity  impact. 

An  increase  of  the  support  size  of  the  projectile  inner  particles  by  a  factor 
of  1.5  has  been  numerically  verified  as  the  optimal  adjustment.  This 
support  size  adjustment  was  done  by  using  a  MATLAB  program  that 
modifies  the  input  file  input_dila.DAT. 


Figure  55.  Projectile  deformation  without  support  size  adjustment. 
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Figure  56.  Projectile  discretization  and  cross  section. 


•  Optimal  smoothing  zone: 

The  optimal  smoothing  zone  for  SNNI  is  defined  by  half  of  the  nodal 
distance,  as  shown  in  Figure  57.  This  parameter  is  implemented  directly  in 
the  latest  version  of  NMAP. 


Figure  57.  Smoothing  zone  for  SNNI 
(smoothing  zone  shaded). 
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•  Summary  of  the  parametric  study: 

Table  8  provides  the  summary  of  guidelines  for  modeling  the  V&V 
penetration  problems  using  NMAP. 


Table  8.  Summary  of  selected  parameters. 


Panel  Discretization 

1.2xl0'3  -1.5*10-3  m  (corresponding  to 
non-dimensional  wave  number  0.045~0.056a) 

Bullet  Discretization 

Nodal  distance  ratio  between  discretization  in  bullet 
and  panel,  1 :  1.4-1.5 

Time-step  Size 

AT  <  1/10*CFL 

Boundary  Layer  Correction 

Move  the  discretized  bullet  surface  inward  by  half  of 
nodal  spacing 

Dynamic  Increase  Factor 
(applied  to  the  bullet) 

DIF  =  1.5  for  Vimpact  >  610  m/sec 

DIF  =  1.2  for  Vimpact  <  610  m/sec 

Support  Size  Adjustment 

Increase  by  a  factor  of  1.5  for  bullet’s  inner  particles 

Smoothing  Zone  Size 

Coefficient  *  Nodal  distance, 

Coefficient  =  0.5 

a  The  non-dimensional  wave  number  is  normalized  by  the  major  wavelength  of  the  system. 
The  non-dimensional  wave  number  smaller  than  0.1  gives  an  accurate  result  when  RKPM 
with  SCNI  is  used  (You  et  al.  2002),  which  is  used  as  a  reference  for  SNNI.  More 
investigation  on  the  stability  condition  will  be  made  in  the  near  future. 

Summary  of  numerical  simulation  of  V&V  problems 

Eighteen  penetration  experiments  were  conducted  by  ERDC  and  are 
compared  to  NMAP  simulations.  Comparison  of  exit  velocities  and  velocity 
reductions  for  each  test  are  given  in  Table  9.  Figures  58  to  60  provide  a 
comparison  of  the  V&V  experimental  data  and  the  NMAP  numerical 
results  for  each  size  of  projectile,  32-grain,  86-grain  and  129-grain, 
respectively.  Figure  61  (a),  (b)  and  (c)  represent  respectively  the  impact 
face,  exit  face  and  hole  of  the  damaged  panel  of  Run  1.  Each  figure  was 
obtained  by  removing  all  particles  that  have  tensile  damage  greater  than 
0.8,  and  velocity  greater  than  5  m/sec.  Since  the  panel  is  discretized  with 
layers  of  particles  through  the  thickness,  the  layer  that  presents  the 
smallest  crater  is  used  to  specify  the  hole  in  the  panel.  For  analysis,  the 
arrows  drawn  on  each  figure  are  measured,  and  extreme  values  are 
considered  as  minimum  and  maximum  dimensions  of  the  impact  crater, 
exit  crater,  and  hole.  These  are  summarized  in  Table  10.  Projectile  and 
panel  mass  loss  and  crater  debris  weights  are  shown  in  Table  11.  The  mass 
loss  criterion  in  the 
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Table  9.  Exit  velocity  comparison. 


Pre-Test 

Post-T  est 

Run# 

Test# 

Panel 

Projectile 

Actual 

Panel 

Thickn 

ess 

(in.) 

Impact 

Velocit 

y 

(ft/s) 

Exit 

Velocity 

(ft/s) 

Velocity 
Reducti 
on  (%) 

Exit  Velocity 
(ft/s) 

NUMERICAL 

Velocity 

Reduction 

(%) 

NUMERICAL 

4 

090806-02 

I 

129  gr 
1/2" 
sphere 

1  5/8 

5636 

2597.00 

53.92% 

2119.42 

62.39% 

8 

090819-01 

1  9/16 

1337 

0.00 

100.00% 

4.56 

99.66% 

15 

090527-06 

86  gr 

7/16" 

sphere 

1  7/16 

2616 

50.00 

98.09% 

-14.27 

100.00% 

trial  1 

091013-03 

1  9/16 

5604 

2246.00 

59.92% 

1873.69 

66.57% 

7 

091014-02 

32  gr 
5/16" 
sphere 

1  9/16 

6718 

177.00 

97.37% 

697.83 

89.61% 

3 

090806-01 

1  9/16 

1985 

0.00 

100.00% 

2.33 

99.88% 

14 

090917-01 

Nominal  1"  Thickness 

129  gr 
1/2" 
sphere 

1 

1582 

75.00 

95.26% 

21.52 

98.64% 

13 

090904-02 

86  gr 

7/16" 

sphere 

1 

3657 

1787.00 

51.13% 

1051.84 

71.24% 

9 

090819-02 

1 

1767 

0.00 

100.00% 

5.51 

99.69% 

12 

090904-01 

1 

884 

0.00 

100.00% 

13.06 

98.52% 

trial  2 

090911-03 

1  1/16 

1892 

0.00 

100.00% 

-7.19 

100.00% 

16 

091014-01 

32  gr 
5/16" 
sphere 

1  1/16 

6521 

1284.00 

80.31% 

2464.57 

62.21% 

11 

090903-01 

1  1/16 

3015 

80.00 

97.35% 

36.22 

98.80% 

6 

090807-02 

Nominal  1/2"  Thickness 

129  gr 
1/2" 
sphere 

9/16 

1516 

654.00 

56.86% 

383.20 

74.72% 

2 

090805-02 

5/8 

377 

0.00 

100.00% 

26.57 

92.95% 

10 

090902-03 

86  gr 
7/16" 
sphere 

9/16 

859 

0.00 

100.00% 

67.03 

92.20% 

1 

090805-01 

32  gr 
5/16" 
sphere 

11/16 

2868 

976.00 

65.97% 

296.59 

89.66% 

5 

090807-01 

17/32 

785 

0.00 

100.00% 

12.20 

98.45% 

Convert  inches  to  millimeters  as  ( •  )x25.4. 


Convert  feet  per  second  to  meters/second  as  ( •  )x0.3048. 
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Figure  58.  Comparison  of  velocity  reduction  for  V&V  tests  using  32-grain  projectiles. 
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Figure  59.  Comparison  of  velocity  reduction  for  V&V  tests  using  86-grain  projectiles. 
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Figure  60.  Comparison  of  velocity  reduction  for  V&V  tests  using  129-grain  projectiles. 


(a)  impact  face  (b)  exit  face  (c)  hole 

Figure  61.  Illustration  of  measurement  of  crater  sizes  and  hole  dimensions. 


panel  is  defined  by  removing  all  particles  that  have  tensile  damage  greater 
than  0.8,  and  velocity  greater  than  5  m/sec.  No  damage  criterion  is 
defined  for  the  projectile.  Therefore,  the  mass  loss  of  the  projectile  is 
measured  by  removing  the  particles  that  separate  from  the  projectile 
during  the  impact  process. 

Velocity  histories  and  damage  patterns  are  presented  for  four  V&V  tests 
(Runs  6,  7, 10  and  16).  Run  6  describes  the  impact  of  a  129-grain  projectile 
of  low  initial  velocity  onto  a  thin  panel  of  nominal  12.7-mm  thickness.  Run  7 
describes  the  impact  of  a  32-grain  projectile  of  high  initial  velocity  onto  a 


Table  10.  Crater  and  hole  diameter  comparison. 
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Table  11.  Projectile  and  panel  mass  loss  and  crater  debris  weight  comparison. 
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thick  panel  of  nominal  38.1-mm  thickness.  Run  10  describes  the  impact  of 
an  86-grain  projectile  of  low  impact  velocity  onto  a  thin  panel  of  nominal 
12.7-mm  thickness.  Run  16  describes  the  impact  of  a  32-grain  projectile  of 
high  impact  velocity  onto  a  panel  of  nominal  25.4-mm  thickness.  The 
velocity  history  curves  of  Runs  6,  7, 10  and  16  are  shown  in  Figures  62,  64, 
66,  and  68,  respectively. 

The  damage  patterns  on  the  impact  face  and  exit  face  of  Runs  6,  7, 10  and 
16  are  also  provided  in  Figures  63,  65,  67,  and  69,  respectively.  Provided  is 
a  comparison  of  the  craters  between  the  numerical  simulation  and  the 
experimental  observation. 

Comparisons  between  macroscale  and  multiscale  AFC  models 

In  this  study,  two  sets  of  V&V  tests  are  adopted  for  comparison.  Since  in 
this  study  the  same  dimensions  for  both  projectiles  and  panels  are  used 
but  with  different  impact  velocities,  the  comparison  is  identified  based  on 
the  impact  velocity  of  the  projectiles. 


Figure  62.  Velocity  history  of  Run  6. 
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Figure  63.  Comparison  of  damage  pattern  for  Run  6. 
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Figure  64.  Velocity  history  of  Run  7. 
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Figure  65.  Comparison  of  damage  pattern  for  Run  7. 


Figure  66.  Velocity  history  of  Run  10. 
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(a)  I Face  -  Experiment 
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Figure  67.  Comparison  of  damage  pattern  for  Run  10. 
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Figure  68.  Velocity  history  of  Run  16. 
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Figure  69.  Comparison  of  damage  pattern  for  Run  16. 


Different  damage  patterns  of  multiscale  and  macroscale  AFC  models  are 
given  in  Figures  70  through  74.  Table  12  lists  the  corresponding  exit 
velocities.  For  the  multiscale  AFC  model,  the  damage  patterns  include 
shear  and  tensile  damage,  while  only  shear  damage  is  considered  in 
macroscale  AFC  model.  Therefore,  the  figures  only  exhibit  the  shear 
damage  pattern  for  both  the  macroscale  and  multiscale  AFC  models. 
Hence,  the  multiscale  AFC  model  is  expected  to  generate  a  larger  damage 
region  than  the  macroscale  AFC  model.  As  the  impact  wave  of  the 
projectile  propagates  to  the  back  of  the  concrete  panel,  then  bounces  back, 
it  causes  tensile  damage  and  therefore  induces  a  larger  damage  zone  in  the 
exit  face.  As  can  be  seen  in  Figures  70  through  74,  the  multiscale  AFC 
model  predicts  larger  craters  than  those  predicted  by  the  macroscale  AFC 
model.  For  Run  7,  the  diagram  of  the  multiscale  AFC  model  seems  to 
generate  a  smaller  crater  than  the  macroscale  one.  This  is  because  the 
tensile  damage  patterns  are  not  drawn  in  the  comparison. 
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(c )  Impac  t  fa  ce  (Mul  d  -  sc  ale)  (d)  Exi  t  face  (Mul  d  -  s  cale) 


Figure  70.  Damage  pattern  of  Run  12. 


(a)  Inpact  face  (Macro-scale)  (b)  Exit  face  (Macro-scale) 


(c)  Impact  face  (Mul d -scale)  (d)  Exit  face  (Multi -scale) 


Figure  71.  Damage  pattern  of  Run  9. 
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(a)  Impact  facs  (Macro-scale) 


(c)  Impact  files  (Multi -seals) 


(b)  Exit  files  (Macro-scale) 


(d)  Exit  face  (Multi -seals) 


Figure  72.  Damage  pattern  of  Run  13. 


(a)  Impact  face  (Macro-scale) 


(c )  Impac  t  fa  ce  (Mul  ti  -  sc  al  e)  (d)  Exi  t  face  (Multi  -s  cal  e) 

Figure  73.  Damage  pattern  of  Run  3. 


(b)  Exit  face  (Macro-scale) 
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(a)  Impact  face  (Macro-scale)  (b)  Exit  face  (M acro-scale) 


(c)  Impact  face  (Multi -sc ale)  (d)  Exit  face  (Multi -scale) 


Figure  74.  Damage  pattern  of  Run  7. 


Table  12.  Comparisons  of  exit  velocity  between  macroscale,  multiscale 
AFC  models  and  experimental  measurement. 


Macroscale 

AFC  Model 
(m/s) 

Multlscale 

AFC  Model 
(m/s) 

Experiments 

(m/s) 

Run  12 

-2.62 

3.98 

0.00 

Run  9 

-0.29 

1.68 

0.00 

Run  13 

350.60 

320.60 

544.68 

Run  3 

0.65 

0.71 

0.00 

Run  7 

129.98 

212.70 

53.95 

Empirical  formula  extracted  from  penetration  V&V  numerical  and 
experimental  data 

Figure  75  illustrates  the  typical  velocity  reduction  characteristics  obtained 
from  the  numerical  simulations  using  an  86-grain  projectile  with  various 
panel  thicknesses  and  impact  velocities.  Here,  a  general  formula  is 
constructed  for  describing  the  projectile  velocity  reduction  based  on 
numerical  and  experimental  data  of  the  V&V  problems.  To  formulate  the 
relationship  among  the  three  parameters,  i.e.,  velocity  and  size  of  the 
projectile  and  the  panel  thickness  two  types  of  criteria  are  considered  for 
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86  gr  bullet 


Figure  75.  Velocity  reduction  rate  with  various  thicknesses  and  impact 

velocities. 


representing  the  projectile  energy.  These  are  velocity  and  kinetic  energy 
criteria,  which  are  proposed  to  quantify  the  behavior  of  the  projectiles. 

•  Velocity  criteria: 

Based  on  the  V&V  problems,  two  observations  are  made  from  the 
numerical  simulations  and  experimental  results. 

o  If  the  impact  velocity  of  the  projectile  is  less  than  a  threshold,  the 
projectile  will  stop  inside  the  panel  or  bounce  back.  In  order  for  the 
projectile  to  perforate  the  panel,  a  minimum  projectile  velocity  is 
required. 

o  Thicker  panels  yield  greater  velocity  reduction,  and  therefore  the 
minimum  velocity  required  to  perforate  is  a  function  of  panel 
thickness. 

Figures  76  to  78  show  the  comparison  between  numerical  and  experimental 
results  using  different  projectile  sizes.  Exit  velocity  versus  normalized 
impact  velocity  (normalized  by  the  panel  thickness)  is  plotted.  Table  13  lists 
the  minimum  normalized  impact  velocity  required  for  perforation,  as 
determined  from  the  numerical  and  experimental  results.  For  the  tests  with 
a  32-grain  projectile,  larger  discrepancy  exists  between  the  numerical  and 
experimental  results,  while  the  simulation  results  for  86-grain  and 
129-grain  projectile  cases  agree  well  with  the  experimental  data. 
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Impact  Velocities  vs  Exit  Velocities 
VeXit(ft/s)  (32-grain  projectile) 

3000  - 


y  =  0.000172  x2  -  0.7797  x  +  766 


Vjmpact/thickness(ft/s/in) 


Figure  76.  Normalized  impact  velocity  vs.  exit  velocity  (32-grain  projectile). 


Figure  77.  Normalized  impact  velocity  vs.  exit  velocity  (86-grain  projectile). 


ERDC/GSL  TR-11-35 


82 


Table  13.  Minimum  required  normalized  velocity. 


(  V  \ 

impact 

32-grain 

86-grain 

129-grain 

thickness  .  w 

V  /  required 

projectile 

projectile 

projectile 

RKPM 

3093  (ft/s/in.) 

1636  (ft/s/in.) 

2049  (ft/s/in.) 

Experiments 

2451(ft/s/in.) 

1678  (ft/s/in.) 

1976  (ft/s/in.) 

Convert  ft/s/in  to  m/s/mm  as  ( •  )x0.012 


•  Kinetic  energy  criterion: 

The  following  assumptions  based  on  the  numerical  and  experimental  data 
are  made  for  the  construction  of  an  empirical  formula  using  kinetic  energy 
criterion. 

o  Energy  loss  occurs  for  the  projectile  during  penetration  of  the 

concrete  panel.  If  the  projectile’s  impact  kinetic  energy  is  less  than  a 
threshold,  it  will  not  perforate  through  the  panel.  Consequently,  a 
minimum  kinetic  energy  is  required  for  perforation, 
o  The  minimum  required  kinetic  energy  is  also  dependent  on  the 
panel  thickness.  The  thicker  the  panel,  the  larger  the  reduction  of 
kinetic  energy. 
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Based  on  the  previous  information,  an  empirical  formula  is  constructed  in 
the  following  form 


Kexit  =  Rd^fL~Kcr  (109) 

where  Kimpact  and  Kmt  are  the  impact  and  exit  kinetic  energy  of  the 
projectile,  Kcr  is  the  critical  kinetic  energy  for  the  projectile  to  perforate 
the  panel,  t  is  the  thickness  of  the  panel  and  Rd  is  the  kinetic  energy 
reduction  factor,  which  is 


Rd=j  (110) 

where  c  is  a  constant  that  may  be  related  to  the  material  properties  of  the 
projectile  and  panel. 

By  least-squares  fitting  of  the  previous  formula  to  the  numerical  and 
experimental  results  as  plotted  in 

Figure  79,  it  can  be  seen  that  the  curves  generated  from  the  numerical  and 
experiment  data  show  similar  trends,  although  discrepancies  exist.  Table 
14  lists  the  values  of  the  characterized  coefficient  c  and  the  critical  kinetic 
energy  Kcritical .  These  procedures  will  be  performed  more  rigorously  based 

on  regression  analysis  or  other  well  developed  mathematical  tools  in  the 
future. 

Summary  and  discussions 

Several  numerical  simulations  were  first  carried  out  to  determine  certain 
optimum  model  parameters.  All  other  V&V  numerical  results  were  based 
on  use  of  these  optimal  values.  The  DIF  is  adopted  since  the  J2  plasticity 
model  for  the  projectile  does  not  consider  rate  effects  to  instantaneously 
adjust  yield  stress  with  the  changing  of  strain  rate.  For  instance,  on  Runs  7 
and  16,  the  experimental  mass  loss  of  the  projectiles  is  about  60%,  while 
numerically  only  a  small  amount  of  mass  is  lost.  Both  cases  have  extremely 
high  impact  velocity  but  with  the  smallest  projectile  size.  Physically,  the 
mass  loss  may  be  due  to  the  melting  of  the  projectile  and  subsequent 
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Impact  Kinetic  Energy  vs  Exit  Kinetic  Energy 

Kexft  (lb-ft) 

2000  -  - 


Figure  79.  Normalized  impact  kinetic  energy  vs.  exit  kinetic  energy. 


Table  14.  Fitted  coefficient  c  and  critical  value  Kcrltlcal . 


c 

is 

A  critical 

RKPM 

0.2778  (in.) 

1321  (lb-ft) 

Experiments 

0.4511  (in.) 

1710  (lb-ft) 

Convert  in  to  mm  as  ( •  )x25.4;  convert  lb-ft  to  J  as  ( •  )xl.3558. 


erosion,  while  the  J2  plasticity  model  does  not  include  the  thermal  effects 
of  metallic  material.  As  the  temperature  of  the  projectile  increases,  the 
strength  of  the  material  should  decline.  This  suggests  that  the  magnitude 
of  DIF  needs  to  be  modified  to  compensate  for  the  thermal  effects.  In  the 
current  NMAP  code,  the  magnitude  of  DIF  employed  is  only  determined 
by  the  impact  velocity.  For  the  impact  velocity  larger  than  6io  m/s,  DIF  of 
1.5  is  used,  while  DIF  of  1.2  is  selected  for  the  impact  velocity  less  than 
610  m/s.  The  larger  DIF  reduces  energy  dissipation  during  perforation 
yielding  an  increase  of  exit  velocity.  Furthermore,  the  damage  model  in  the 
projectile  is  not  considered  in  the  present  study.  These  are  the  reasons  why 
the  numerical  results  exhibit  higher  exit  velocity  than  the  experimental 
ones  for  these  two  cases.  If  the  thermal  and  damage  effects  are  considered 
in  the  numerical  models,  more  dissipated  energy  will  be  encountered  in 
the  numerical  simulation,  and  lower  numerical  exit  velocities  are  expected. 
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Empirical  formulas  are  proposed  to  quantify  the  velocity  reduction  of 
projectile.  As  can  be  seen  from  Figures  76  to  78,  the  numerical  V&V  results 
show  a  good  agreement  with  the  experiments.  Based  on  the  velocity 
criterion,  the  numerical  results  show  similar  trends  in  comparison  with 
the  experimental  data  except  for  the  cases  with  32-grain  projectile.  As 
aforementioned,  this  can  be  attributed  to  the  fact  that  thermal  effects  and 
projectile  damage  mechanisms  are  not  considered  in  the  simulation 
models. 

For  the  material  model  of  concrete,  a  multiscale  AFC  model  is  proposed. 
For  the  proposed  multiscale  AFC  model,  tensile  and  shear  damage  effects 
are  considered  separately,  while  in  the  original  macroscale  AFC  model, 
only  shear  damage  effect  is  included  in  the  numerical  algorithm.  As  can  be 
seen  from  the  damage  patterns,  the  multiscale  AFC  model  exhibits  larger 
damage  zones  than  the  macroscale  model  due  to  the  consideration  of 
tensile  damage. 
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6  Conclusions 

A  multiscale  RKPM  formulation  for  impact  and  penetration  modeling  was 
developed,  with  specific  focus  on  methods  to  accurately  model  key 
phenomena  such  as  multi-body  contact,  fast-evolving  material  damage, 
large  material  deformation,  and  large  material  flow.  An  enhanced  semi- 
Lagrangian  formulation  was  introduced  to  address  the  breakdown  of  the 
Lagrangian  approximation  when  conditions  such  as  free  surface  formation 
or  closure  occur,  which  commonly  result  from  material  fracture  during 
penetration  events.  Temporal  stability  of  the  semi-Lagrangian  formulation 
was  also  analyzed,  and  a  strong  dependence  was  shown  between  stability 
criteria  and  the  velocity  gradient.  Dependence  of  the  stability  criteria  on 
velocity  gradient  is  of  particular  importance  in  impact  problems,  where 
very  large  gradients  are  encountered  in  the  impact  region.  For  strain 
calculation,  a  stabilized  non-conforming  nodal  integration  scheme  was 
introduced.  The  SNNI  scheme  was  derived  from  the  similar  stabilized 
conforming  nodal  integration  method,  but  with  simpler  non-conforming 
strain  smoothing  domains  that  do  not  require  continuous  reconstruction 
in  the  semi-Lagrangian  formulation. 

To  model  multi-body  contact,  an  enhanced  kernel  contact  algorithm  was 
developed.  Contact  conditions  were  detected  based  on  kernel  interaction 
between  nodes  of  different  bodies,  allowing  body  interaction  to  be 
captured  naturally  within  the  RKPM  framework.  As  such,  a  priori 
definition  of  contact  surfaces  was  not  required,  which  is  a  significant 
enhancement  for  problems  where  extensive  material  separation  and  debris 
formation  occur.  A  technique  using  the  level  set  method  was  introduced  to 
determine  surface  normal  at  a  contact  surface.  The  level  set  method 
provides  a  more  accurate  determination  of  the  surface  normal  as 
compared  to  using  a  point-to-point  technique.  Decomposition  of  the 
internal  force  between  two  nodes  was  used  to  identify  the  normal  force 
acting  between  them,  and  tangential  forces  were  subsequently  modified 
using  the  Coulomb  friction  law.  Under  separation  conditions,  internal 
forces  were  modified  to  prevent  non-physical  tensile  forces  that  were 
shown  to  give  erroneous  pull-back  conditions. 

A  new  multiscale  formulation  for  damage  evolution  was  developed,  where 
continuum-scale  damage  evolution  was  linked  to  microstructure  failure. 
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First  a  microcell  analysis  was  introduced  to  explicitly  model  crack 
propagation  in  the  material  microstructure.  Relationships  between  stress, 
strain  and  Helmholtz  free  energy  of  the  microcell  and  homogenized 
continuum  were  then  used  to  derive  a  new  continuum-scale  damage 
evolution  function,  which  was  termed  the  microcrack-informed  damage 
model.  It  was  shown  that  this  multiscale  approach  could  be  used  to  derive 
microcrack-informed  damage  laws  ranging  from  simple  one-  and  two- 
parameter  models  to  full-tensorial  damage  descriptions.  The  numerical  size 
effect  between  the  microcell  and  structural-scale  models  was  also  studied, 
and  a  scaling  law  was  presented  to  obtain  mesh  objectivity.  Notably,  this 
multiscale  formulation  provides  a  framework  in  which  numerical 
experiments  at  the  microscale  are  used  to  determine  continuum-scale 
parameters  that  were  otherwise  based  on  phenomenological  laws.  As  such, 
much  more  accurate  descriptions  of  these  parameters  (such  as  the  damage 
evolution  function)  are  provided,  yielding  a  desired  capability  in  terms  of 
structural-scale  analysis  and  multiscale  material  design. 

New  formulations  developed  in  the  project  were  implemented  into  the  3-D 
NMAP  parallel  code.  The  microcrack  informed  damage  model  was 
implemented  into  the  ERDC-developed  Advanced  Fundamental  Concrete 
model,  and  was  used  to  define  tensile  damage  evolution  in  the  CorTuf 
ultra  high-strength  concrete  material.  This  MIDM-enhanced  version  of  the 
AFC  model  was  also  implemented  into  NMAP  for  the  analysis  of  projectile 
penetration  into  ultra  high-strength  concrete  targets. 

The  new  formulations  and  NMAP  code  implementation  were  evaluated 
through  18  V&V  tests  that  were  based  on  penetration  experiments 
conducted  by  ERDC.  A  small  set  of  numerical  experiments  were  initially 
performed  to  determine  optimal  values  for  certain  model  parameters.  Key 
findings  from  the  parameter  optimization  study  included  discretization 
guidelines,  critical  time-step  criteria,  dynamic  increase  factor  criteria,  and 
recommendations  for  boundary  layer  corrections,  projectile  support  size 
adjustment  and  smoothing  zone  geometry.  In  terms  of  the  V&V  test 
results,  generally  good  agreement  was  found  from  comparison  of  the 
numerical  results  and  experimental  data.  With  regard  to  percent  velocity 
reduction  of  the  projectile,  error  of  10%  or  less  was  observed  in  14  of  the 
tests.  Error  in  the  remaining  four  tests  was  less  than  24%,  where  these 
four  cases  coincided  with  strong  overmatch  conditions  in  the  experiments. 
Several  of  the  V&V  tests  (Runs  3,  7,  9, 12  and  13)  were  selected  to  compare 
results  using  the  original  AFC  model  and  the  MIDM-enhanced  version. 
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Each  case  was  run  using  both  AFC  model  versions,  and  comparisons  were 
made  in  terms  of  exit  velocity  and  panel  damage.  It  was  found  that  the 
MIDM-enhanced  AFC  model  generally  predicted  higher  exit  velocity, 
larger  craters  and  more  damage  than  the  original  AFC  model.  These 
results  were  in-keeping  with  a  more  accurate  tensile  failure  description 
provided  by  the  multiscale  formulation.  Lastly,  numerical  and 
experimental  results  were  used  to  develop  empirical  equations  for 
predicting  terminal  ballistic  conditions.  Velocity  criteria  or  kinetic  energy 
criteria  were  used  as  a  basis  for  the  empirical  expressions,  and  either  exit 
velocity  or  exit  kinetic  energy  was  estimated  as  a  function  of  impact 
velocity  or  energy  normalized  by  the  panel  thickness.  These  empirical 
estimators  were  shown  to  provide  reasonable  predictions  over  the  given 
range  of  impact  conditions. 
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Appendix  A:  Derivation  of  Influence  Tensor 
Using  Asymptotic  Expansion  Based 
Homogenization 

An  alternative  approach  for  relating  microscopic  and  macroscopic 
variables  is  the  asymptotic  expansion  based  homogenization  method. 
Consider  the  model  problem  described  by  coarse-scale  coordinate  x  and 
fine-scale  coordinate  y  as  shown  in 

Figure  14.  The  two  length  scales  are  related  with  a  small  parameter,  s,  as 

1 

y  =  — x  (Al) 

e 

The  spatial  derivatives  of  a  function  O'  with  superscript  s  denoting  the 
combined  (total)  coarse-fine  scale  characters  is  expressed  as 

VxOs(x)  =  VxO(x,y)  =  VxO(x,y) + - V  yO(x,y)  (A2) 

£ 

Based  on  asymptotic  expansion,  the  total  displacement  vector  is  expanded 
as 


u'(x)  =  J>iuI'1(x,y)  (A3) 

i=0 

Substituting  (A3)  into  Equation  41  and  considering  (At)  and  (A2), 
expansion  of  the  strain  tensor  is  obtained  as 

00 

es  =  Vxu£  =  £kelk]  (A4) 

k=- 1 


where 


=  VyU[0] 


eLfcJ(u)  =  Vxu[fc]  +  V£u[fc+1] 


(A5) 


k>  0 


(A6) 
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Introducing  the  elasticity  tensor,  the  stress  expansion  is 

oo 

o'  =  C* :  e'  =  £  eVa 

k=- 1 


where 


C  :  e 


[*] 


(A7) 


(A8) 


Substituting  expansion  (A8)  into  the  equilibrium  condition  Equation  36, 
the  leading  order  equilibrium  equations  is  obtained 


Vy  -a[~1]  =0 

(A9) 

Vx •  a[k]  +Vv •  o[/c+1]  =  0  k>- 1 

x  y 

(A10) 

The  formal  solutions  of  these  equations  were  discussed  in  the  literature 
(Cheng  1992;  Fish  et  al.  1997;  Guedes  and  Kikuchi  et  al.  1990)  and 
outlined  are  only  the  important  results  here.  The  first  and  second  order 
formal  solutions  for  displacement  are  expressed  as 

u[0]=v[0](x)  (All) 


u[1]=v[1](x)  +  x(y):Vsxv[0](x) 


(A12) 


where  v[0](x)  and  v[1](x)  are  the  coarse  and  fine  scale  solution  functions  of 
u,  and  x(y)  is  the  third  order  characteristic  tensor  function  of  the 
microscopic  cell  (Bakhvalov  and  Panasenko  1989). 

Consider  the  truncation  of  the  strain  and  stress  expansions  to  two  scales 


e£  =-e[  1] +e[0]  + . »-e[  1]+e[0]  (A13) 

£  £ 


o£  =— o[  1] +a[0]  + . a[  1]  +o[0]  (A14) 

£  £ 


where: 
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e[-1]  = Vyii[0] = Vy  v[0]  (x)  =  0  ( A15) 

o[_1] =C£ :  e[_1] =C£ :  v;u[0]  =C£ :  V^v[0]  (x)  =  0  (A16) 

e[0]  =  Vxu[0]  +Vyu[1]  =  [l  +  Vyx(y)]:  Vxv[0](x)  (A17) 

a[0] =C£ :  e[0]  =  C£ :  [l  +  V*x(y)] :  Vsx v[0]  (x)  (A18) 

Substituting  (Ai5~i8)  into  (Ai3~i4),  results  in  the  following 

e£  «  e[0]  =  Vxu[0]  +  Vyu[1] = [l  +  V*x(y)] :  Vsxu[0]  (x)  (A19) 

a£  w  o[0]  =  C£ :  [l  +  Vsyx(y)] :  V£u[0]  (x)  (A20) 


where  I  denotes  the  fourth-order  identity  tensor.  Substituting  stress 
expression  (A20)  into  equilibrium  Equation  36  and  the  cohesive  traction 
Equation  39,  the  following  equations  result  for  solving  the  third-order 

tensor,  x(y) 

vy-(c£:v;x(y))  =  0  in  ny  (A21) 

[c£ :  Vyx(y):  Vxv[0](x)]-n  =  -[c£ :  Vxv[0](x)]  n  +  h  on  Yc  (A22) 

Comparing  Equation  52  and  (A19),  the  influence  tensor  can  be  expressive 
by  the  characteristic  tensor  as 


A£=-v;x(y) 


(A23) 
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Appendix  B:  Derivation  of  Energy  Bridging 
Equation  Using  Asymptotic  Expansion  Based 
Homogenization 

Here  we  show  that  the  energy  bridging  Equation  69  can  be  obtained  by  an 
asymptotic  expansion  based  method.  Substituting  Equations  (A19)  and 
(A20)  into  the  microscopic  free  energy  defined  in  Equation  64  yields 

We  =  ia* .  e*  =  ie* .  Ce .  Qe  =  i(vsxu[0]  +  v;u[1]) :  C£ :  (v*u[0]  +  v;u[1]) 

2  2  2  (B1 

=  — Vxu[0] :  CE :  Vxu[0]  +-Vs„u[1] :  Ce :  V*u[1]  +  V^u[0] :  C£ :  Vsu[1] 


Integrating  over  the  microscopic  cell,  the  following  is  obtained 


f  wedn  =  1 

f— Vxu[0] : 

:Ce  :Vxu [0]}d£l 

2 

+Xjv;uIU 

\ 

:  C£ :  V£u[1] 

/ 

dn  +  (Vxu[0] :  C£ :  Vsyu[1])dQ 

(B2) 


The  three  terms  on  the  right  hand  side  of  Equation  (B2)  are  rearranged  as 
follows.  The  first  term  is  expressed  as 


-  V£u[0] :  C£ 
12 


V;u 


[0] 


dtl  =  ^vsxu 
2  x 


[0] 


(c£):Vxu[0] 


(B3) 


The  second  term  can  be  further  manipulated  by  considering  k  =  - 1  in 
Equation  (A10)  to  yield 

V-o[“1]+V-o[0]  =  0  inO  (B4) 

Further  considering  the  crack  surface  traction  condition  as 

o[0]  n  =  h  on  Tc  (B5) 

and  substituting  Equation  (A16)  in  Equation  (B4),  the  following  results 
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Vy  -o[0]  =  0  (B6) 

Multiplying  Equation  (B6)  by  u[1]  and  integrating  it  over  Q  yields 

f  u[1]-v-o[0]dn  =  0  (B7) 

y 

By  integration  by  parts,  yields  the  following 

f  u[1]-V-o[0]dQ=  f  o[0]  •  V  •  u[1]dn  -  f  u[1]  •  o[0]  •  ndS  =  0  (B8) 

J  ny  y  J  aY  y  J  rc 

where  Yc  is  the  crack  surface  within  the  microscopic  cell.  Considering  the 
boundary  condition  in  Equation  (B5)  and  the  symmetry  of  stress  tensor, 
o[0] ,  the  following  is  shown 

f  o[0] :  Vju[  =  -  f  u[1]-h dS  (B9) 

JnY  y  Jrc 

Substituting  Equation  (A18)  into  Equation  (B9),  one  obtains 

Vsxu[0] :  C£ :  Vyu[11df2  +  V^u[1] :  C£ :  Vsyu[11d£2  =  u[1]  •  h dS  (BIO) 

The  second  term  on  the  right  side  of  Equation  (B2)  is  then  obtained.  For 
the  third  term  on  the  right  hand  side  of  Equation  (B2),  considering 
Equation  (A12) 

f  Vxu[0] :  Ce :  vyi]dft  =  Vsxu[0] :  [  f  Ce :  Vyx(y)rfnl :  Vxu[0]  (Bll) 

Jily  y  [Jfly  ^  J 

Substituting  Equation  (A20)  into  the  homogenized  stress  expressed  in 
Equation  49  yields 

d  =  (o£)  =  (C£ :  fl ■ +  V*X(y)l) :  V£u[0](x) 

\  /  \  l  J/  (B12) 

=  (C£) :  Vxu[0](x)  +  (C£ :  Vsyx(y)) :  V£u[0](x) 

Comparing  Equation  (B11)  with  (B12)  gives 
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Then  substituting  Equations  (B3),  (Bio),  and  (B13)  into  Equation  (B2),  the 
following  is  obtained 

— (  f  x¥ed£l  +  —  f  u[1]  hds|  =  — d:  V*u[0]  (B14) 

Vy  (4  Hy  2  4  r c  J  2 


It  is  observed  that  V'u[0]  is  the  macroscopic  strain,  e ,  that  is  imposed  on 

the  microscopic  cell  by  Equation  67.  By  substituting  Equation  60  in 
Equation  (B14),  the  following  is  finally  obtained 

W  =  —  f  WedQ  +  -  f  u[1]-hcis]  (B15) 

2drc  J 
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Appendix  C:  Two-Variable  Damage  Model 

According  to  the  continuum  damage  theory,  the  general  expression  of  the 
fourth-order  damage  representation  is 

a  =  (l-D):C0:ee  (Cl) 

where  C0  is  the  initial  elastic  stiffness  tensor,  D  is  the  fourth  order 
damage  tensor,  and  ze  =e  —  sp  and  zp  is  the  plastic  strain.  The  effective 
stress  is  defined  as 


C0 :  £e 


(C2) 


Hence,  the  damage  model  in  Equation  (Cl)  can  be  expressed  as 


o  =  (l-D):o0 

(C3) 

To  account  for  the  unilateral  effect,  the  positive-negative  decomposition  of 
the  effective  stress  tensor  is  defined  as 

O0  =0+  +o0 

(C4) 

o 

D 

+ 

CL 

II 

+  o 

D 

(C5) 

o 

D 

i 

CL 

II 

+  o 

D 

1 

o 

D 

II 

+  o 

D 

(C6) 

where  the  fourth-order  projection  tensors  P+  and  P 

are  (Faria  et  al.  1998) 

P+  =  p.  ®p;  <8>pf  0P; 

i 

(C7) 

P  =I-P+ 

(C8) 

in  which  I  is  the  fourth-order  identity,  d;  and  p,  are  the  ith  eigenvalue 
and  the  corresponding  eigenvector  for  the  effective  stress  tensor  o0 ,  and 
H(x)  is  the  Heaviside  function 
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H(x)  = 


x>0 

x<0 


(C9) 


Correspondingly,  two  damage  scalars,  d+  and  d~ ,  are  introduced  to 
describe  the  damage  of  materials  under  tension  and  compression 
respectively.  According  to  the  principles  of  thermodynamics,  the  state  of 
an  ensemble  could  be  described  by  using  the  definition  of  Helmholtz  free 
energy  (HFE),  which  is  expressed  by  using  the  state  variables  and  internal 
variables.  The  total  elastoplastic  HFE  is  defined  as 

ip  =  ip(ee,K,d+,d~)  (CIO) 


where  k  is  the  plastic  variables.  Decomposing  the  total  HFE  into  the 
elastic  and  the  plastic  components,  yields 

yj(se ,  k,  d+ ,  d )  =  ipe  (ee ,  d+ ,  d )  +  ipp  (ee ,  k,  d+ ,  d )  (Cll) 


Neglecting  the  plastic  strain  under  tension,  the  following  plastic  HFE  is 
described 


xpptze,K,d+,d )  =  ipptee,K,d -)  =  (1  -d~)vp0  (C12) 

The  elastic  HFE  is  further  decomposed  as 

y/(e  e,d+,d-)  =  ipe+(ze,d+)  +  ipe+(se,d-) 

=a-d+w0+(sn+a-d-w0-(ee)  1 


Here  the  subscripts  “e”  and  “p”  refer  to  “elastic”  and  “plastic”  components, 
respectively,  and  the  subscript  “o”  refers  to  the  “initial”  state. 

According  to  the  second  principle  of  thermodynamics,  the  following 
equation  is  obtained 

0  =  Itt-- (1  -  d+K  +  (1  -  d-  )0-  -  (I  -  d+P+  -  d-p-) :  o0  (C14) 

oz 

The  damage  energy  release  rate  (DERR)  can  then  be  obtained  as 

Y±  _  dty_ 

dd± 


(C15) 
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